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Abstract 
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In  this  project,  the  components  of  a  Stirling  space  nuclear  power  system  were  modeled 
using  MATLAB®  and  Simulink®  under  both  steady  state  and  transient  conditions.  Using 
information  provided  through  NASA’s  Glenn  Research  Center,  the  Department  of 
Energy’s  Naval  Reactors,  and  literature  from  previous  work  in  nuclear  space  power 
systems,  the  system  architecture  was  developed.  At  steady  state,  the  space  power  system 
was  designed  to  produce  200  kW(e)  and  consisted  of  a  NaK  cooled  fast  reactor  linked  to 
four  Stirling  power  converters  via  two  NaK  heat  transfer  loops  connected  by  heat 
exchangers.  A  third  NaK  heat  transfer  loop  was  used  to  transport  the  waste  heat  of  the 
Stirling  converters  to  radiator  panels,  which  radiate  the  heat  into  ambient  space.  Heat 
pipes  were  modeled  which  linked  the  secondary  heat  transfer  loop  to  the  Stirling 
converters  and  the  third  heat  transfer  loop  to  the  radiator  panels.  An  alternator  was  used 
to  convert  the  thermal  power  input  of  the  Stirling  converters  into  electrical  power.  The 
model  was  based  on  combining  the  point  reactor  kinetic  equations  with  transient  energy, 
force,  volume,  and  mass  balances  for  system  components  outside  the  reactor.  The  model 
was  tested  with  dynamic  perturbations  that  included  small  reactivity  changes  in  the 
reactor,  Stirling  converter  failures,  effects  of  planetary  albedo,  and  electrical  resistance 
load  changes.  At  steady  state,  the  model  produced  the  expected  temperatures,  pressures, 
flow  rates,  and  heat  flows.  The  dynamic  simulations  indicated  that  the  overall  system 
would  be  load  following  and  stable  to  external  perturbations.  However,  in  cases  that 
varied  the  electrical  load  resistance,  the  thermal  efficiency  of  the  Stirling  converters  were 


significantly  affected. 
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pw  wick  porosity 

s  radiator  emissivity 

sr  regenerator  effectiveness 

r|  heat  exchanger  efficiency 

r|ait  alternator  efficiency 

T]conv  instantaneous  converter  efficiency 
X  precursor  decay  constant 

p  fluid  viscosity 

pm  coolant  viscosity 

p  reactivity 

p  fluid  density 

Pf  fuel  density 

pgr  graphite  density 

pm  coolant  density 

pn  flow  density  in  space  n 

PNaK  NaK  coolant  density 

p0  initial  reactivity 

cr  Stephan-Boltzman  constant 

t  time  delay  for  the  coolant  to  travel  from  the  reactor  to  the 

heat  exchanger 
co  angular  frequency 

C,  damping  ratio 
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"In  conclusion  then,  the  navigation  of  interplanetary  space  depends  for  its  solution  on  the  problem  of 
atomic  disintegration...  Thus,  something  impossible  will  probably  be  accomplished  through  something  else 
which  has  always  been  held  equally  impossible,  but  which  remains  so  no  longer. " 

-  Robert  H.  Goddard,  October  3,  1907 


1.  Background 

While  many  components  come  together  to  make  a  spacecraft,  the  component  that  has 
the  most  significant  impact  on  the  spacecraft  is  its  power  system.  If  a  spacecraft  was 
compared  to  a  human,  its  power  system  would  be  described  as  its  cardiovascular  system. 
All  other  functions  of  the  spacecraft  rely  on  it  and  the  spacecraft  would  be  inoperable 
without  it.  This  is  best  shown  by  the  recent  Japanese  spacecraft  ADEOS-2,  Advanced 
Earth  Observing  Satellite,  which  has  been  deemed  “all  but  lost”  after  its  power  system 
suffered  a  catastrophic  failure.1  Therefore,  space  power  systems  must  not  only  be 
designed  to  facilitate  the  operations  of  the  spacecraft,  but  they  must  also  be  designed  to 
withstand  the  variety  of  stimulus  that  the  space  environment  has  to  offer.  Large  changes 
in  current  or  voltage  within  the  power  system  components  can  lead  to  decreased 
performance  or  failure  of  the  subsystems.  The  failure  or  degradation  of  these 
components  can  render  the  spacecraft  incapable  of  performing  its  mission. 

In  the  recent  past,  the  space  community  has  consistently  used  photovoltaic  and 
chemical  means  of  providing  power.  These  power  systems  have  performed  very  well  for 
their  low  power  missions  near  Earth.  However,  their  application  to  missions  farther  from 
the  Earth  or  missions  requiring  higher  power  levels  is  very  limited.  Current  state  of  the 
art  solar  cells,  such  as  Spectrolab’s  Ultra  Triple  Junction  solar  cells,  can  reach 
efficiencies  of  almost  thirty  percent.  '  With  a  solar  flux  at  the  Earth  of  approximately 
1360  W/m2,  a  solar  panel  that  is  directly  facing  the  sun  at  all  times  could  theoretically 


produce  power  densities  of  around  400  W/m2.  To  have  a  solar  array  capable  of 


Solar  Flux  as  a  Function  of  Distance  from  the  Sun 


13 


Figure  1-1.  Solar  energy  flux  relative  to  the  Earth.  (AU=astronautical  unit,  the  distance  between  the 

centers  of  the  Earth  and  Sun,  1.495x10s  km) 

producing  200  kW(e)  of  power,  which  is  the  approximate  power  level  of  the  Prometheus 
spacecraft,4  a  solar  array  of  500  square  meters  would  be  needed.  An  array  of  this  size 
would  be  the  equivalent  size  of  a  tenth  of  a  football  field.  This  calculation  neglects  the 
degradation  of  the  panels  over  time,  the  fraction  of  light  reflected  off  the  panels,  the 
possibility  of  the  Earth  or  other  objects  blocking  the  sunlight,  and  the  extra  panels  needed 
for  a  power  margin.  The  energy  flux  from  the  sun  is  also  inversely  proportional  to  the 
square  of  the  distance  from  the  Sun.  So  for  a  spacecraft  orbiting  Jupiter,  which  is  five 
times  farther  from  the  Sun  than  the  Earth,  twenty-five  times  more  solar  panel  area  is 
needed  to  produce  the  same  amount  of  power  (Figure  1-1). 5  This  means  that  a  solar 
panel  area  equal  to  that  of  about  three  footballs  fields  would  be  required.  Improvements 
in  efficiency  and  mass  alone  will  not  solve  these  problems  since  the  huge  array  areas 


required  for  high  power  will  create  obstacles  in  structural  integrity,  stowing  for  launch, 
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deployment  in  orbit,  and  sun  pointing  that  are  far  from  being  solved  with  current 
technology.6  The  effects  of  the  Sun  also  prevent  the  use  of  solar  power  for  missions 
closer  to  the  sun.  The  increase  in  energy  and  radiation  from  the  Sun  destroy  solar  panel 
materials  when  spacecraft  get  too  close. 


Figure  1-2.  Space  power  sources  as  a  function  of  mission  power  level  and  duration.  [8] 

In  an  effort  to  increase  the  capabilities  of  space  exploration,  the  space  community  has 

turned  to  nuclear  power  as  a  source  of  continuous  power  for  deep  space  missions.  For 
lower  power  missions,  radioisotope  thermal  generators  (RTG’s)  have  been  extremely 
successful.  RTG’s  use  the  thermal  energy  released  by  radioactive  decay  to  produce 
electricity.  Missions  that  have  used  RTG’s  include  Pioneer,  Voyager,  Viking,  and  the 
more  recent  Cassini  spacecraft.  The  Pioneer  10  spacecraft  was  still  running  on  its  RTG 
after  30  years  in  space.7  The  load  produced  by  various  space  power  sources  versus 

o 

mission  duration  is  given  in  Figure  1-2.  As  can  be  seen,  RTG’s  are  limited  to  around  2 
kW(e)  if  they  are  to  be  cost  efficient. 


In  order  to  accomplish  deep  space  missions  that  require  higher  power  levels,  the  space 
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community  has  turned  to  fission  reactors  for  the  solution.  This  is  because  fission  of  one 
kilogram  of  uranium-235  produces  500,000  times  the  amount  of  energy  released  by  the 
decay  of  one  kilogram  Plutonium-238  (used  in  RTG’s)  over  ten  years.9  Starting  as  early 
as  the  mid- 1 950’ s,  the  United  States  explored  fission  reactors  as  a  space  power  source 
with  the  Systems  for  Nuclear  Auxiliary  Power  (SNAP)  program.  This  program  resulted 
in  the  only  reactor  flown  in  space  by  the  United  States.  Also,  Russia  has  been  using 
small  fission  reactors  as  space  power  sources  since  the  early  1960’s.10  In  the  late  1970’s, 
the  United  Nations  Working  Group  on  the  Use  of  Nuclear  Power  Sources  in  Outer  Space 
recognized  the  advantages  of  nuclear  power  in  space  and  recommended  its  use  for 


Figure  1-3.  SP-100  surface  power  concept  (left)  and  nuclear  electric  propulsion  spacecraft  concept  (right).  [12] 
technical  pursuits.11  Shortly  thereafter,  the  United  States  reinitiated  its  pursuit  of  nuclear 

fission  power  in  space  with  the  SP-100  program  (Figure  1-3).  The  SP-100  program 

lasted  until  the  early  1990’s  and  consisted  of  studies  for  various  nuclear  fission  power 

spacecraft  configurations  and  missions.  The  studies  provided  positive  insight  into  the 

use  of  nuclear  fission  power  in  space  and  recommended  several  concepts.  Since  there 

was  never  a  mission  requirement  for  the  SP-100  program,  the  program  was  canceled 

1  T 

without  ever  completing  a  space  mission. 


Today,  even  with  the  advances  in  solar  panel  materials  and  efficiency  over  the  past 
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decade,  NASA  has  come  to  the  conclusion  that  nuclear  fission  power  is  necessary  for 
high  power,  deep  space  missions.  NASA  recently  initiated  the  Prometheus  Project  to 


Figure  1-4.  The  conceptual  design  of  the  JIMO 
nuclear  electric  propulsion  spacecraft.  [14] 

make  the  vision  of  space  nuclear  power  and  propulsion  into  a  reality  by  developing 
missions  along  with  the  technology.  The  first  proposed  mission  for  the  project  was  the 
Jupiter  Icy  Moons  Orbiter  (JIMO,  Figure  1-4),  which  would  use  a  nuclear  powered, 
electric  propulsion  spacecraft  to  investigate  three  of  Jupiter’s  moons  in  great  detail.14 
NASA  has  recently  reconsidered  JIMO  and  postponed  the  mission  due  to  the  high  risks 
involved.15  Instead,  the  Prometheus  Program  will  develop  a  lower  risk  mission  for  the 
initial  Prometheus  undertaking.  NASA  is  pursuing  this  lower  risk  and  lower  cost 
program  to  increase  the  likelihood  of  success.  One  of  the  proposals  for  this  initial 
mission  includes  a  surface  power  system  for  the  Moon  or  Mars,  similar  to  the  SP-100 
program. 

With  the  technology  provided  by  NASA,  a  mission  to  Mars  can  be  achieved  three 
times  faster  than  with  current  chemical  propellants.16  For  propulsion,  chemical 


propellants  are  the  standard  fuel  for  launch  into  Earth  orbits.  While  this  will  most  likely 
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continue  to  be  the  case  for  future  missions,  propulsion  from  Earth  orbits  to  distant 
locations  of  the  solar  system  become  costly  or  time  consuming  using  chemical 
propellants.  For  example,  the  mass  of  the  current  Atlas  V  500  series  rocket  with  Centaur 
upper  stage  is  91.5%  propellant.  This  leaves  little  room  for  improvement  in  terms  of 
launching  a  spacecraft  farther  into  the  solar  system.  Either  the  mass  of  the  spacecraft  has 
to  be  reduced  by  removing  components  or  more  fuel  has  to  be  added.  If  the  spacecraft’s 
time  of  flight  is  to  be  reduced,  even  more  fuel  must  be  added.  This  is  due  to  the  fact  that 
specific  impulse,  which  measures  how  efficiently  a  rocket  can  achieve  a  velocity  change 
for  a  given  mass  of  propellant,  is  four  to  twelve  times  less  for  chemical  propellants  than 
for  the  electrostatic  ion  engines  proposed  for  the  JIMO  mission.  As  a  result,  spacecraft 
using  chemical  propulsion  are  limited  to  having  only  10-20%  of  its  mass  as  the  payload, 
while  nuclear  power  and  propulsion  are  able  to  raise  this  number  to  forty  percent.19  With 

nuclear  power,  ion  engines  (Figure  1-5)  would  have  the  capability  of  continuous 

20 

operation  without  prohibiting  the  use  of  other  instruments  or  systems  on  the  spacecraft. 
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Figure  1-5.  Conceptual  diagram  of  electric  propulsion.  [20] 
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However,  the  nuclear  reactor  is  just  one  of  several  system  components  that  will  make 
up  the  space  power  system.  Another  major  system  component  is  the  energy  converter, 
the  object  which  will  convert  the  thermal  energy  of  the  nuclear  reactor  into  electrical 
power.  Most  highly  developed  space  nuclear  power  system  studies  have  involved 
thermoelectric  converters,  which  use  the  electrical  potential  created  by  the  temperature 
gradient  across  a  material  to  produce  electricity.  However,  these  converters  have  very 
low  efficiency  (~  5%)  and  NASA  is  currently  considering  the  mechanical  alternatives 
using  power  converters  based  on  Brayton  and  Stirling  cycles. 

In  a  NASA  contractor  report,  the  performance  of  Stirling  and  Brayton  converters  in 
the  SP-100  power  system  is  compared.““  The  results  of  the  report  assert  that  when  an  SP- 
100  reactor  is  coupled  with  Brayton  or  Stirling  cycle  engines  with  the  same  power 
requirement,  the  Stirling  cycle  option  has  a  greater  net  cycle  thermal  efficiency  and  lasts 
longer  (Table  1-1).  This  means  that  Stirling  cycle  engines  have  the  potential  for  longer 


Parameter 

Brayton 

Stirling 

Net  Power  (kW) 

550 

550 

Gross  Power  (kW) 

582 

596 

Reactor  Thermal  Power  (kW) 

2,309 

1,850 

Net  Cycle  Thermal  Efficiency  (%) 

23.8 

29.7 

Possible  Carnot  Efficiency  (%)* 

68.77 

52.41 

Full  Power  Lifetime  (yrs) 

7.6 

9.6 

Operating/Redundant  Converters 

3/1 

3/1 

Temperatures  (K) 

Reactor  Outlet 

1,355 

1,355 

Average  Maximum 

1,300 

1,265 

Average  Minimum 

406 

602 

Radiator  Area  (m2) 

590 

210 

System  Mass  (kg),  1000  Vdc 

15,768 

14,932 

System  Mass  (kg),  1000  Vrms,  1  kHz 

14,903 

15,217 

Table  1-1.  SP-100  Brayton  vs.  Stirling  comparison.  [22] 
*Field  added  by  the  author. 


mission  life  or  lower  launch  weight.  Additionally,  the  Brayton  cycle  option  has  received 
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significant  attention  in  terms  of  analysis  and  modeling.  Thus  a  study  of  the  Stirling 
option  could  produce  some  meaningful  results  in  an  area  that  has  received  less  attention. 
This  is  especially  true  for  a  model  of  large  Stirling  engines  such  as  those  used  for  the  SP- 
100  design.  Creating  a  model  using  this  size  of  a  Stirling  engine  and  coupling  the  engine 
to  a  nuclear  reactor  are  both  very  unique  undertakings.  For  this  project,  a  computer 
model  of  a  Stirling  nuclear  power  system  was  created  using  MATLAB®  and  Simulink® 
software.  The  model  development  completed  the  following  objectives: 

(1)  Developed  time-variant  computer  models  for  the  thermodynamic, 
mechanical,  and  electrical  processes  of  the  reactor,  Stirling  engines, 
radiators,  alternators,  and  heat  transfer  components. 

(2)  Used  empirical  data  from  Naval  Reactors  and  NASA  to  determine  system 
parameters  (sizes,  ratios,  etc.). 

(3)  Validated  the  model  mathematically  and  physically  (engineering). 

(4)  Tested  the  system’s  stability  by  modeling  the  response  to  possible 
perturbations  (i.e.  engine  failure,  reactivity  changes). 

(5)  Increased  system  performance  characteristics  such  as  efficiency  or  specific 
power. 

The  most  significant  objectives  were  (1)  and  (2).  All  other  objectives  were  secondary 
and  were  based  on  the  thorough  completion  of  the  first  two  objectives. 

The  system  studied  consists  of  several  components  which  have  been  previously 
investigated  in-depth.  While  the  behavior  of  the  individual  components  of  this  system  is 
well  studied  and  understood,  the  resulting  behavior  of  the  combination  of  the  non-linear 
systems  has  never  been  undertaken  before  now.  Since  Stirling  converters  are  being 
considered  as  one  of  the  main  options  for  future  space  power  system  development,  this 
study  provides  knowledge  for  the  possible  development  of  these  future  systems. 


2.  System  Overview 

The  process  of  nuclear  electric  power  (NEP)  can  be  simply  stated  as  the  conversion  of 
thermal  energy  released  by  nuclear  reactions  to  electrical  energy.  For  this  research,  the 
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Figure  2-1.  General  NEP  structure.  [20] 

three  sections  of  the  NEP  system  that  were  modeled  and  evaluated  are  the  reactor,  power 
conversion,  and  heat  rejection  components.  This  section  gives  an  overview  of  the 
components  that  are  included  within  these  sections  for  a  Stirling  nuclear  space  power 
system. 

The  reactor  chosen  for  this  research  is  a  liquid  metal  fast  reactor  (LMFR).  A  LMFR 
consists  of  an  original  fuel  source  of  enriched  uranium,  which  is  uranium  that  has  a 
higher  than  naturally  occurring  percentage  of  uranium-235.  This  enriched  uranium  also 
contains  uranium-238,  which  can  absorb  fast  neutrons  to  produce  plutonium-239.  In 
unconstrained  fission,  a  neutron  is  absorbed  by  a  fissile  nucleus.  The  nucleus  becomes 
unstable  and  splits,  releasing  energy,  neutrons,  and  fission  products  (Figure2-2)2  For 
uranium-235  and  plutonium-239,  2.6  and  2.98  neutrons  are  released,  respectively,  on 
average  in  a  fast  reactor.25 


A  fast  reactor  is  a  system  in  which  neutrons  maintain  an  energy  near  the  magnitude 
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they  are  produced.  To  accomplish  this,  materials  which  contain  low  atomic  mass 


An  example  of  one  of  the  many 
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Figure  2-2.  Fission  of  uranium-235.  In  a  fast  reactor,  the  average  number  of  neutrons  released 

increases  to  2.6.  [24] 

elements,  such  as  water,  must  not  be  present.  Liquid  metals  work  well  since  they  have 


relatively  high  atomic  masses.  For  this  research,  sodium-potassium  metal  (NaK)  is  used 


as  the  liquid  metal  coolant.  Sodium-potassium  allows  the  reactor  to  achieve  high  power 


densities,  defined  as  power  produced  per  unit  volume,  due  to  its  high  conductivity 
combined  with  its  relatively  high  boiling  point.” 


The  heat  transfer  from  the  reactor  to  the  Stirling  engines  is  facilitated  by  two  sodium- 


potassium  heat  transfer  loops  as  shown  in  Figure  2-3.  Two  fluid  loops  are  used  to  avoid 


material  degradation  issues  that  may  result  from  the  radioactive  isotopes  created  in  the 


primary  loop  when  coolant  passes  through  the  reactor.  For  example,  sodium  will  absorb 
neutrons  to  form  the  beta-gamma  emitter  Na.  By  using  two  loops  linked  with  a  heat 
exchanger,  the  radioactive  isotopes  can  be  contained  within  the  primary  loop  and  the 
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4  Stirling  Converters _ _ _ 

Secondary  NaK  Loop  \\  U 
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Figure  2-3.  Schematic  of  the  Stirling  space  nuclear  power  system  (not  to  scale), 
reactor  while  a  non-radioactive  secondary  loop  can  transfer  heat  directly  to  the  Stirling 

converters.  This  two  loop  configuration  is  commonly  used  with  a  LMFR. 

The  problem  of  moving  the  heat  from  the  second  heat  transfer  loop  to  the  Stirling 
engine  was  investigated  during  the  design  process  for  the  SP-100  power  system.  A 
method  chosen  by  the  SP-100  program  is  shown  in  Figure  2-4  and  has  been  adopted  for 
this  project.26  This  method  involves  using  heat  pipes  with  evaporator  ends  exposed  to  the 
secondary  loop  coolant  flow  and  condenser  ends  integrated  into  the  geometry  of  the 
Stirling  hot  end  heat  exchanger.  A  heat  pipe  is  a  device  that  uses  a  working  fluid 
saturated  between  liquid  and  gas  phases  to  transfer  heat.  A  saturated  fluid  is  a  fluid  that 


Figure  2-4.  Design  to  transfer  heat  from  the  liquid  metal  of  the  second  heat  transfer  loop  to  the 

heater  of  the  Stirling  engine.26 


is  in  the  process  of  changing  phases,  such  as  liquid  to  gas.  While  changing  phases,  the 


23 


working  fluid  maintains  an  almost  constant  temperature  as  the  energy  it  receives  is  used 
to  overcome  the  latent  heat  of  vaporization  and  complete  the  phase  change.  The  working 
fluid  absorbs  energy  at  one  end  of  the  heat  pipe  and  releases  energy  at  the  other,  which 
keeps  the  working  fluid  temperature  within  the  saturated  region  (Figure  2-5).  The 
constant  temperature  of  the  saturated  fluid  results  in  a  small  temperature  drop  across  the 
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Figure  2-5.  Basic  process  of  a  heat  pipe,  left,  and  a  cross-section  of  an  actual  heat  pipe,  right.  [27,  46] 
heat  pipe. 

A  Stirling  converter  is  a  device  capable  of  using  the  temperature  difference  provided 
by  a  heat  source  and  a  heat  sink  to  produce  work  via  a  gaseous  working  fluid.  The  ideal 
Stirling  cycle  consists  of  four  processes:  isothermal  compression,  isometric  heat  addition, 
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Figure  2-6.  Ideal  Stirling  cycle  with  single  piston  and  displacer. 
(www.grc.nasa.gOv/VVVVW/tmsh/stirling/Slidepage4.html) 


isothermal  expansion,  and  isometric  heat  rejection  (Figure  2-6).  This  ideal  cycle  has  the 
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capability  of  producing  Carnot  efficiency,  which  is  the  maximum  theoretical  efficiency  a 
thermodynamic  cycle  can  have  for  a  given  hot  to  cold  fluid  temperature  ratio  (Th/Tc). 
However,  for  the  practical  application  of  Stirling  power  converters,  this  ideal  cycle  is  not 
beneficial  and  can  never  be  completely  achieved.  In  order  to  model  such  an  ideal  cycle 
where  net  work  is  maximized,  the  ideal  Stirling  P-V  relationship  shown  in  Figure  2-7 
requires  the  frequency  of  the  cycle  to  approach  zero.  Yet,  the  purpose  of  the  Stirling 
converter  is  to  produce  power,  which  is  the  product  of  the  net  work  created  per  cycle  and 
the  frequency.  If  the  frequency  of  the  converter  approaches  zero,  the  power  of  the 
converter  approaches  zero  and  converter  power  is  sacrificed  for  efficiency.  So,  in  reality, 
the  Stirling  converter  is  operated  at  some  frequency  and  the  pressure  volume  diagram 
becomes  less  ideal  in  shape  (Figure  2-7). 29,30 

The  main  components  of  the  Stirling  converter  are  the  heater,  regenerator,  cooler, 
displacer,  power  piston,  and  alternator.  The  heater  and  cooler  provide  a  continuous  heat 
source  and  sink,  respectively,  for  the  Stirling  converter.  The  regenerator  adds  efficiency 
to  the  cycle  by  cooling  the  hot  gas  as  it  moves  from  the  expansion  space  to  the 
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Volume  (cc) 


Figure  2-7.  Comparison  of  an  ideal  (left)  and  a  practical  (right)  pressure-volume  diagram  for  a 
Stirling  converter.  The  net  work  of  the  cycle  is  the  area  within  the  curve. 


compression  space  and  by  heating  the  cold  gas  as  it  moves  from  the  compression  space  to 
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the  expansion  space.  The  displacer  is  a  moving  component  that  separates  the  cold 
compression  space  from  the  hot  expansion  space.  It  is  the  pressure  and  temperature 


Figure  2-8.  Conceptual  drawing  of  JIMO  propulsion,  radar,  and  communication  systems  [31]. 
gradients  across  the  displacer  that  create  force  imbalances  on  the  moving  components  of 

the  Stirling  converters,  the  piston  and  displacer.  These  forces  convert  the  thermal  energy 

of  the  working  fluid  into  the  motion  of  the  piston  or  displacer.  The  mechanical  energy,  or 

motion,  of  the  piston  is  then  converted  into  electrical  energy  using  a  linear  alternator. 

The  linear  alternator  consists  of  magnets  that  are  linked  to  the  piston  and  coiled  wires  that 

are  linked  to  the  Stirling  casing.  The  motion  of  the  magnets  produces  a  voltage  within 

the  coils  based  on  the  velocity  of  the  piston.  This  voltage  then  creates  current  and 

electric  power  based  on  the  impedance  of  the  associated  electrical  circuit.  The  resulting 

electrical  power  can  be  used  for  scientific  instruments,  energy  storage,  or  any  other 

device  that  requires  electricity.  For  the  JIMO  mission,  this  would  include  the  Prometheus 

spacecraft’s  radar  and  electric  propulsion  units,  and  other  spacecraft  subsystems  (Figure 

2-8). 31 

In  order  for  Stirling  power  converters  to  convert  thermal  energy  into  electrical  power, 
there  must  be  a  heat  source  and  a  heat  sink.  The  heat  source  for  the  Stirling  converters  is 


the  heat  generated  by  the  nuclear  reactor.  The  heat  sink  for  the  Stirling  converters  is  deep 
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space,  and  the  heat  rejection  is  accomplished  using  heat  pipes  in  combination  with 
radiator  panels.  The  excess  heat  of  the  Stirling  cycle  is  transferred  to  a  liquid  water  loop 
by  forced  convection.  Heat  pipes  transfer  the  heat  from  the  fluid  loop  to  the  radiators. 
Due  to  the  small  temperature  drop  across  heat  pipes,  they  provide  an  equalized  energy 
distribution  to  the  radiators  that  helps  maximize  the  radiation  of  heat  into  space  for  a 
given  radiator  geometry,  which  in  turn  minimizes  the  amount  of  radiator  surface  area 
needed.  Since  radiators  are  the  second  largest  source  of  system  mass,  minimizing  their 
mass  is  an  important  factor  in  the  overall  system  design,  but  was  beyond  the  scope  of  this 
project. 


Volume 

Accumulator 


Figure  2-9.  Diagram  of  the  SP-100  general  system  design.  The  system  modeled  by  this  research  is  of 

the  same  general  design.  [22] 

The  following  sections  (3-5)  provide  a  detailed  description  of  the  model  chosen  and 
the  applicable  assumptions  and  parameters  for  each  component  of  the  Stirling  space 
nuclear  power  system.  The  mathematics  and  science  behind  each  model  are  elaborated 
upon  and  justified.  Section  3  discusses  the  reactor  and  the  connected  heat  transfer  loops. 


Section  4  describes  the  Stirling  converters  and  Section  5  details  the  heat  rejection  system. 
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Then,  the  component  and  integrated  component  simulation  results  are  presented  and 
discussed  in  Section  6.  The  final  section.  Section  7,  concludes  with  a  brief  summary  of 
the  work  accomplished,  conclusions  made,  and  a  recommendation  for  future  work. 
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3.  Reactor,  Heat  Exchanger,  and  Heat  Pipe  Models 

3.1  Nuclear  Reactor 

The  reactor  that  is  modeled  for  this  project  is  a  600  kW(t)  fast  fission  reactor 
consisting  of  uranium  nitride  fuel  and  sodium  potassium  coolant.  Its  dynamic  behavior  is 
modeled  using  the  one  group  point  kinetics  equations,  which  treats  the  reactor  as  a  point 
source  of  power  and  ignores  the  loss  of  fuel  mass  over  time.  The  latter  assumption  does 
not  produce  significant  error  since  tests  are  run  over  relatively  short  time  periods  with 
respect  to  fuel  consumption.  The  equations  used  also  simplify  the  many  possible 
combinations  of  fission  products  into  one  average  group.  The  point  kinetics  equations 
use  the  properties  of  the  reactor’s  fuel,  coolant,  and  geometry  to  model  the  reactor’s 
power,  component  temperatures,  and  other  factors  over  time. 


clP_ 

dt 


— — —P  +  XPE 


dPE 

dt 


-P-XPE 

A 


3-1 


3-2 


P  =  Po+af(Tf-Tfo)  +  am  -  Tmo  )  3'3 

where:  P  =  reactor  power 

PE  =  reactor  precursor  power  equivalence 
p  =  reactivity  (p0  =  initial  reactivity) 

[)  =  delayed  neutron  fraction 
A  =  neutron  generation  time 
X  =  precursor  decay  constant 
a/  =  fuel  temperature  coefficient 
am  =  coolant  temperature  coefficient 
Tm  =  coolant  temperature 
Tf  =  average  fuel  core  temperature 
Tmo  =  equilibrium  coolant  temperature 
Tfo  =  equilibrium  fuel  temperature 


Equations  3-1  and  3-2  of  the  point  kinetics  equations  relate  the  rate  of  change  of  the 
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reactor’s  power  and  precursor  power  equivalence  to  the  time  dependent  variables  of 
reactivity  ( p ),  power  (P),  and  precursor  power  equivalence  (PE)  and  the  reactor  constants 
of  delayed  neutron  fraction,  neutron  generation  time,  and  precursor  decay  constant.  The 
power  of  the  reactor  is  a  result  of  fission  and  the  precursor  power  equivalence  is  a 
function  of  fission  product  decay.  The  reactivity  of  a  reactor  is  a  function  of  the  neutron 
multiplication  factor,  k,  where  k  is  the  ratio  of  the  number  of  neutrons  in  the  present 
generation  and  the  number  of  neutrons  in  the  preceding  generation  (Equations  3-4  and  3- 
5).  For  a  reactor  producing  a  steady  power  level,  the  number  of  neutrons  in  one 
generation  produces  an  equal  amount  of  neutrons  in  the  next  generation.  This  state 
corresponds  to  a  multiplication  factor  value  of  one  and  a  reactivity  value  of  zero. 


P 


k- 1 
k 


#  of  neutrons  in  present  generation 
#  of  neutrons  in  the  preceding  generation 


3-4 


3-5 


Equation  3-3  determines  the  reactor’s  reactivity  at  a  given  time  according  to  the  initial 
values  of  the  reactivity,  fuel  temperature,  coolant  temperature  and  the  values  of  fuel  and 
coolant  temperature  and  temperature  coefficient.  The  temperature  coefficients  of  the  fuel 
and  coolant  use  a  linear  approximation  to  characterize  the  change  in  the  reactor’s 
reactivity  due  to  temperature  change.  Values  for  the  temperature  coefficients  are  a 
function  of  the  materials  that  make  up  the  fuel  and  coolant.  The  temperature  coefficient 
plays  a  large  role  in  the  stability  of  the  reactor  since  it  defines  the  response  of  the  reactor 
due  to  temperature  change.  The  temperatures  of  the  fuel  and  the  coolant  affect  their 
density  and  nuclear  cross  sections,  which  in  turn  affect  the  overall  reactivity  of  the 


reactor.  To  have  a  self-stabilizing  reactor,  a  negative  fuel  temperature  coefficient  is 
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required.  A  negative  fuel  temperature  coefficient  decreases  the  reactivity  as  temperature 
rises  and  increases  the  reactivity  as  temperature  decreases.  Therefore,  a  negative  fuel 
temperature  coefficient  results  in  the  reactor  naturally  trying  to  maintain  its  current  power 
level.  On  the  other  hand,  a  positive  fuel  temperature  coefficient  will  cause  the  reactor  to 
diverge  from  its  current  power  level. 

The  neutron  generation  time  is  typically  about  the  same  value  for  all  fast  reactors. 
The  delayed  neutron  fraction  and  decay  constant  are  properties  of  Uranium-235.  Using 
these  values  and  the  initial  condition  of  zero  reactivity,  the  steady  state  reactor  precursor 
power  equivalence  can  be  determined  for  a  600  kW(t)  reactor  power  (Equations  3-1  and 
3-2).  This  power  level  was  calculated  using  the  200  kW(e)  power  requirement  and  the 
assumption  that  the  Stirling  converters  have  about  thirty-three  percent  efficiency.  The 
power  level  of  the  reactor  changes  the  temperatures  of  the  fuel  and  coolant  according  to 
relationships  that  are  dependent  on  the  reactors  geometry,  materials,  and  coolant  flow 
rate.  For  a  square  lattice  circular  fuel  pin  geometry  (Figure  3-1),  these  relationships  can 


Figure  3-1.  Circular  fuel  pin  geometry  with  square  lattice. 

be  expressed  by  Equations  3-6  through  3-9. 34  The  state  of  the  different  system 
components  is  monitored  by  calculating  temperature  as  a  function  of  time.  This  method 
makes  it  easier  to  discover  model  errors  and  predict  component  interaction. 
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where:  pf  =  fuel  density 

pm  =  coolant  density 
Cf  =  fuel  specific  heat 
Cpm  =  coolant  specific  heat 

Cp,ni  =  coolant  specific  heat  at  reactor  inlet  temperature 

Cpm2  =  coolant  specific  heat  at  reactor  outlet  temperature 

N  =  number  of  fuel  rod 

Lf  =  length  of  a  fuel  rod 

kf  =  fuel  thermal  conductivity 

rf  =  radius  of  fuel  pin 

h  =  convection  coefficient 

mm  =  coolant  mass  flow  rate 

Tf  =  average  fuel  core  temperature 

Tm  =  (Tmi  +  Tw2)/2  =  average  coolant  temperature  in  reactor 
Tml  =  coolant  temperature  of  reactor  inlet 
Tm2  =  coolant  temperature  at  reactor  outlet 
Vf  =  fuel  volume  in  reactor 
Vm  =  coolant  volume  in  reactor 

Originally,  a  lumped  capacitance  model  was  considered,  which  neglects  the  thermal 
temperature  gradients  within  the  fuel  pin  and  assumes  that  the  fuel  is  one  uniform 
temperature.  However,  for  the  lumped  capacitance  method  to  be  valid,  the  Biot  number 

ic 

has  to  be  much  less  than  one  (BicO.l).  The  Biot  number  is  calculated  using  Equations 
3-8  and  3-9. 


Bi 
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—  Nu 

de 
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3-9 


where: 


k,„  =  conductivity  of  the  coolant 
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de  = 


4  A 


equivalent  diameter  of  the  flow  passage 


A  =  flow  area  of  the  unit  cell  (Figure  3-2) 

Pw  =  wetted  perimeter  (2jirf  for  this  geometry) 

Nu  =  Nusselt  number 

Since  the  Biot  number  was  close  to  one  hundred,  the  thermal  resistance  of  the  fuel  rod 
was  included  by  using  the  geometry  shown  in  Figure  3-1.  The  temperature  of  the  fuel  pin 
surface  can  be  eliminated  by  mathematical  manipulation  of  the  conduction  and 
convection  equations,  resulting  in  the  form  of  Equations  3-6  and  3-7.  The  Nusselt 
number  in  Equation  3-10  is  valid  when  the  condition  in  Equation  3-11  is  met.  The  high 
conductivity  of  liquid  metals  results  in  their  ability  to  assume  a  constant  Nusselt  number 
based  on  the  pitch  to  diameter  ratio  of  the  fuel  lattice  if  the  Peclet  number  is  below  the 
critical  Peclet  number.  The  critical  Peclet  number  for  this  problem  was  460,  which  is 
significantly  greater  than  NaK’s  Peclet  number  of  17  in  this  case. 


Nu  =  7.0  +  4.24(p/  /  D/)1'52 

3-10 

if 

Pe  <  Pecr 

3-11 

where:  py  =  fuel  pin  pitch  (Figure  3-2) 

Df=  2iy  =  fuel  pin  diameter  (Figure  3-2) 

Pe  =  RePr  =  Peclet  number 

Re  =  Reynolds  number,  Re  =  pmude  /  pm 

Pr  =  Prandtl  number,  Pr  =  Cpmpm  /  km 

Pecr  =  critical  Peclet  number 
Cpm  =  coolant  specific  heat 
pm  =  coolant  viscosity 

Coolant  enters  the  reactor  at  a  temperature  lower  than  the  fuel’s  temperature.  Energy 
is  then  transferred  from  the  fuel  to  the  coolant  as  the  coolant  flows  along  the  length  of  the 


fuel  rod.  The  coolant  mass  flow  rate  and  temperature  difference  across  the  fuel 


determines  the  thermal  energy  that  the  fuel  adds  to  the  coolant.  This  interaction  can  be 
modeled  using  the  conservation  of  energy  in  the  reactor  coolant  (Equation  3-7). 

The  geometry  of  the  reactor  consists  of  cylindrical  fuel  rods  arranged  in  a  square 
lattice  (Figure3-2  and  3-3).  The  reactor  can  be  modeled  as  a  series  of  unit  cells,  shown  in 
Figure  3-2.  The  number  of  fuel  pins,  fuel  mass,  coolant  mass,  and  steady  state 


Figure  3-2.  (left)  A  fuel  rod  conducts  the  heat  from  the  nuclear  reactions  to  its  surface  by  the 
relationship  Qcond=4reLfkf(TrTs),  where  Ts  is  the  fuel  surface  temperature.  The  heat  is  then  convected  to 
the  coolant  according  to  Qconv=2h7rrfLf(Ts-Tm).  (right)  Unit  cell  showing  fuel  pitch,  Pf,  and  diameter,  D. 


^ out  pm2^ m2 


Figure  3-3.  The  reactor’s  geometry  consists  of  cylindrical  fuel  pins  arranged  in  a  square  lattice.  The 
coolant  enters  from  one  end,  is  heated  by  the  fuel  pins,  and  exits  the  other  end  of  the  reactor. 
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temperature  differences  across  the  reactor  were  based  on  estimates  provided  by  Naval 
Reactors.4  Using  the  100  degrees  Kelvin  coolant  temperature  drop  across  the  reactor 
provided  by  Naval  Reactors,  the  mass  flow  rate  through  the  reactor  was  calculated  with 
Equation  3-12. 


P 

m  = - 


3-12 


^  pm  ,^'m\ ) 

Where:  P  =  reactor  power 

cpm  =  (cpmi+cpm2)/2  =  average  coolant  specific  heat  in  reactor 
Tm2  -  Tmi  =  coolant  temperature  drop  across  the  reactor 

The  size  of  the  fuel  pins  was  then  calculated  by  breaking  the  reactor  down  into  unit  cells 
(Figure3-2).  The  wetted  area  for  heat  transfer  per  unit  cell  was  found  to  be  the  sum  of 
four  quarters  of  a  cylinder,  or  the  equivalent  of  one  cylinder.  Based  upon  the  range  of 
50-70  centimeters  for  the  fuel  pin  length  provided  by  Naval  Reactors,4  Equation  3-13  was 
then  used  to  determine  the  radius  of  the  fuel  pin  for  given  pitch  to  diameter  ratio.  This 
equation  describes  the  relationship  between  the  power  per  fuel  rod  and  the  fuel  pins 
geometry  and  material  properties. 
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4nLfkf 
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The  coolant  convection  heat  transfer  coefficient,  h,  was  calculated  using  Equations  3-9 
and  3-10.  The  value  for  the  pitch  to  diameter  ratio  was  iterated  until  a  fuel  temperature 
was  found  that  was  well  below  the  melting  point  of  the  uranium  nitride  fuel.  The 
equilibrium  reactor  parameters  for  the  model  found  by  this  process  are  listed  in  the  table 


below. 
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Sub-System  Parameter 

Steady  State  Value 

Reactor  Power,  kW(t) 

600 

Reactor  Precursor  Power  Equivalence,  GW 

50.6 

Reactor  Reactivity 

0 

Fuel  Temperature,  K 

1097.6 

Average  Reactor  Coolant  Temperature,  K 

950 

Reactor  Coolant  Inlet  Temperature,  K 

900 

Reactor  Coolant  Outlet  Temperature,  K 

1000 

Table  3-1.  Reactor  Model  Parameters 


3.2  Heat  Transfer  Loops 

After  the  heat  from  the  reactor  is  transferred  to  the  coolant,  the  coolant  must  travel  to  a 
heat  exchanger  via  a  transfer  loop.  Within  the  primary  transfer  loop,  the  higher 
temperature  coolant  passes  through  piping  from  the  reactor  to  the  heat  exchanger.  Upon 
exiting  the  heat  exchanger,  the  cooler  fluid  passes  through  a  pipe  from  the  heat  exchanger 
back  to  the  reactor.  Due  to  the  rate  at  which  the  coolant  is  being  pumped  in  the  loop,  the 
axial  conduction  and  heat  losses  along  the  length  of  the  pipe  were  ignored. 

The  two  main  losses  in  the  transfer  loops  is  the  radiation  of  heat  from  the  pipe  into 
space  and  the  loss  in  pressure  due  to  the  friction  between  the  fluid  and  the  walls  of  the 
pipe.  The  loss  due  to  radiation  has  been  previously  calculated  by  Naval  Reactors  and 
found  to  be  three  to  four  orders  of  magnitude  lower  than  the  power  of  the  reactor.4 
Consequently,  radiation  losses  were  neglected  in  this  study.  The  pressure  loss  of  the  fluid 
in  the  pipes  is  overcome  by  pumps.  The  power  required  to  run  these  pumps  is  a  parasitic 
loss  that  decreases  the  total  available  power  produced  by  the  Stirling  converters.  In  any 
future  system  design,  it  will  have  to  be  decided  whether  to  overcome  this  loss  by 


increasing  the  power  level  of  the  reactor  or  to  decrease  the  power  available  to  run  the 


desired  scientific  equipment  or  other  electrical  loads,  but  this  decision  is  beyond  the 
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scope  of  this  project. 

The  NaK  heat  transfer  loops  were  modeled  with  a  time  delay.  The  temperature  of  the 
coolant  coming  out  of  the  reactor  becomes  the  temperature  of  the  coolant  entering  the 
heat  exchanger  after  the  length  of  time  required  for  the  coolant  to  transverse  the  piping 
between  the  reactor  and  the  heat  exchanger.  This  effect  was  modeled  using  the  following 
equation. 


dt 
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where:  Tbi  =  coolant  temperature  at  heat  exchanger  entrance 

Tm 2  =  coolant  temperature  leaving  the  reactor 
x  =  time  delay  for  the  coolant  to  travel  from  the  reactor  to  the 
heat  exchanger 

Additionally,  Simulink®  provided  a  time  delay  tool  which  modeled  the  time  lag  exactly 
by  holding  values  constant  and  outputting  them  at  a  specified  length  of  time  later.  All  the 
other  sections  of  transfer  loops  were  modeled  in  this  same  manner. 


3.3  Heat  Exchanger  Linking  Primary  and  Secondary  Transfer  Loops 

In  order  to  transfer  the  heat  from  the  primary  transfer  loop  to  the  secondary  transfer 
loop,  a  shell  and  tube  heat  exchanger  as  shown  in  Figure  3-4  was  assumed  as  this  type  of 


secondary  fluid  (NaK) 


primary 

fluid 

(NaK) 


Figure  3-4.  Shell  and  tube  heat  exchanger,  with  tube  pitch,  pt. 


heat  exchanger  is  commonly  employed  in  reactor  systems.  Due  to  the  high  conductivity 
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of  the  primary  coolant  (NaK),  a  one  pass  heat  exchanger  was  assumed  to  be  sufficient. 
The  primary  loop  coolant  enters  the  tubes  of  the  heat  exchanger  at  the  temperature  it  left 
the  reactor.  While  in  the  heat  exchanger,  the  primary  coolant  transfers  heat  to  the 
secondary  coolant  (also  NaK)  that  is  flowing  around  the  outside  of  the  tubes  via 
conduction  through  the  heat  exchanger  tubes.  The  Biot  number  was  calculated  for  a 
single  tube  and  found  to  be  on  the  order  of  one  hundred,  which  means  that  the  heat 
exchanger  tubes  cannot  be  treated  with  a  lumped  capacitance  model.  Therefore,  the  heat 
transfer  rate  between  the  primary  and  secondary  fluid  was  calculated  using  the 
relationships  given  in  Equations  3-15  and  3-16. 
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where:  Q  =  heat  transfer  rate  between  primary  and  secondary  fluids 

UA  =  overall  heat  transfer  coefficient  times  the  equivalent 
surface  area  of  the  heat  exchanger 
AT i m  =  log  mean  temperature  difference 
hp  =  convection  coefficient  for  heat  exchanger  primary  coolant 
hs  =  convection  coefficient  for  exchanger  secondary  coolant 
Ap  =  nd,L,Nt  =  heat  transfer  area  for  primary  coolant 
As  =  nd0LtN,  =  heat  transfer  area  for  secondary  coolant 
di  =  heat  exchanger  tube  inner  diameter 
d„  =  heat  exchanger  tube  outer  diameter 
k,  =  tube  conductivity 
L,  =  tube  length 

N,  =  number  of  tubes  in  heat  exchanger 


Since  the  temperature  of  the  primary  and  secondary  coolant  vary  along  the  length  of  the 


heat  exchanger,  the  log  mean  temperature  difference  was  used  to  determine  the  heat 

re 

transfer  rate.  The  log  mean  temperature  difference  is  defined  in  Equation  3-17. 
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Figure  3-5.  Individual  tube  of  the  heat  exchanger.  The  hotter  primary  coolant  passes 
through  the  inside  of  the  tube  while  the  colder  secondary  coolant  flows  around  the  tube. 
( riim  =  primary  mass  flow  rate,  ms  =  secondary  mass  flow  rate) 
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where:  Th,i  =  hotter  fluid  inlet  temperature 

Th,0  =  hotter  fluid  outlet  temperature 
TcJ  =  colder  fluid  inlet  temperature 
TCj0  =  colder  fluid  outlet  temperature 

The  heat  transfer  coefficients  for  the  primary  and  secondary  fluid  were  assumed  constant 
due  to  the  high  thermal  conductivity  of  the  liquid  metal  coolant.  However,  the  flow 
geometry  is  different  than  that  of  the  reactor,  so  Equations  3-18  and  3-19  were  used  to 
determine  the  heat  transfer  coefficients.36 


Nuj  =4.82+0.697 


Nu2 


7+4.24 


1.52 


o  7 


3-18 

3-19 


where:  Nui  =  Nusselt  number  for  first  heat  exchanger  primary  fluid 

Nu2  =  Nusselt  number  for  first  heat  exchanger  secondary  fluid 
pt  =  heat  exchanger  tube  pitch 


Using  a  process  similar  to  that  described  above  for  the  reactor,  the  dimensions  of  the 


39 


heat  exchanger  were  determined.  However,  the  only  available  value  for  the  heat 
exchanger  was  the  inner  diameter  of  the  tubes,  which  was  provided  by  a  NASA 
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contractor  report  on  the  SP-100  system  and  used  as  a  guideline  for  heat  exchanger  size.' 
At  a  steady  state  condition,  the  temperature  drop  of  the  primary  fluid  in  the  heat 
exchanger  is  equal  to  its  temperature  rise  in  the  reactor.  Since  the  secondary  loop  also 
uses  sodium-potassium,  the  secondary  fluid  also  has  the  same  temperature  drop.  This  left 
the  thickness  of  the  tubes,  pitch  of  the  tubes,  and  length  of  the  tubes  as  variables.  By 
using  the  detailed  sketches  of  the  SP-100  design  provided  by  the  NASA  contractor  report, 
the  length  of  the  heat  exchanger  was  chosen  to  be  a  half  meter.  The  pitch  and  thickness 
of  the  heat  exchanger  tubes  were  then  adjusted  to  provide  a  high  heat  exchanger 
efficiency,  which  was  calculated  using  Equation  3-20. 


T  —T 

c.o  h,o 
7]  = - 

Th,i  ~  Tho 

where:  ;/  =  exchanger  efficiency 

Tco  =  cold  side  exit  temperature 
Tu  o  =  hot  side  exit  temperature 
T i,  j  =  hot  side  inlet  temperature 


3-20 


In  other  words,  the  temperature  difference  between  the  hot  temperature  of  the  primary 


and  secondary  loops  was  minimized.  Using  this  process,  the  efficiency  of  the  exchanger 


was  increased  to  7 1  percent,  but  further  increases  in  efficiency  were  determined  to  be  a 


secondary  goal  compared  to  the  primary  effort  of  completing  the  system  model. 
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3.4  Heat  Pipes  Linking  Secondary  Transfer  Loop  to  the  Stirling 

Engines 

The  method  employed  to  link  the  secondary  heat  transfer  loops  to  the  Stirling 
converters  consisted  of  placing  the  evaporator  end  of  heat  pipes  into  the  liquid  metal  flow 
of  the  secondary  loop  and  the  condenser  ends  into  the  metal  structure.  This  configuration 
was  the  most  documented  method  found  by  the  author.26  The  heat  pipe  exterior  and  wick 
are  composed  of  titanium  and  the  working  fluid  is  potassium.  These  material  choices 
were  based  on  the  recent  work  done  by  the  Jet  Propulsion  Laboratory  in  the  area  of  high 
temperature  heat  pipe  applications.37 


Evaporator  End 


Condenser 

Adiabatic  Center  £ncj 


Flow  In 


Heat  Pipes 


Cross 

Section 


Flow  Out 

Figure  3-6.  Conceptual  sketch  of  the  transfer  loop-Stirling  heat  pipe  interface  (not  to  scale). 


The  heat  pipe  heat  transfer  of  the  heat  pipes  was  calculated  by  determining  the  thermal 

TO 

resistances  of  the  various  portions  of  the  heat  pipe.  The  thermal  resistance  of  the 
working  fluid,  potassium,  was  assumed  to  be  negligible  compared  to  the  resistances  of 
the  pipe  and  the  wick  at  the  condenser  and  evaporator.  The  thermal  resistance  of  the  pipe 

io 

and  wick  of  the  evaporator  and  condenser  are  given  by  Equations  3-21  through  3-24. 


ln(dop/dow) 

27rLekp 


3-21 


where:  Rpe  =  thermal  resistance  of  the  pipe  at  the  evaporator 

dop  =  heat  pipe  outer  diameter 
dow  =  outer  diameter  of  the  wick 


Le  =  length  of  the  evaporator 
kp  =  conductivity  of  the  pipe 
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R, ..  = 


Hdow/dJ 
2i tL  k  „ 

e  eff 


3-22 


where:  Rwe  =  thermal  resistance  of  the  wick  at  the  evaporator 

djW  =  inner  diameter  of  the  wick 
keff  =  effective  conductivity  of  the  wick 


Hdop  / dow) 

2tt  Lckp 


3-23 


where:  Rpc  =  thermal  resistance  of  the  pipe  at  the  condenser 

Lc  =  length  of  the  condenser 


R,  „  = 


Hdow/dJ 
2itL  k  „ 

c  eff 


3-24 


where:  Rwc  =  thermal  resistance  of  the  wick  at  the  condenser 


The  value  for  the  effective  thermal  conductivity  of  the  wick  is  a  combination  of  both 
the  wick  material  and  that  of  the  working  fluid.  The  wicks  of  the  heat  pipes  chosen  for 
this  project  consist  of  a  woven  screen  mesh  of  the  same  material  as  the  pipe  and  their 

io 

effective  thermal  conductivity  is  found  using  the  following  relationship. 

k1((k1+kpHl-/3w)(k1-kp)) 

eff  (k1+kp)+(l-/?w)(krkp) 


where:  ki  =  thermal  conductivity  of  the  working  fluid 

fkw  =  wick  porosity  =  0.65 

The  heat  transfer  from  the  liquid  metal  of  the  secondary  loop  to  the  heater  of  the 
Stirling  converter  was  calculated  by  including  the  thermal  resistance  of  the  heat  transfer 


from  the  liquid  metal  to  the  pipes. 


42 


dQkP 


AT 

ER 


^hp\p^evap 


T<i  Twh _ 

Rpe  +  Rwe  +  Rpc  +  R, 


3-26 


where:  dQi,p  =  heat  transfer  through  the  heat  pipes 

Td  =  average  liquid  metal  temperature  at  the  evaporator 
Twh  =  heater  wall  temperature 
NhP  =  number  of  heat  pipes 

hip  =  heat  transfer  coefficient  from  the  liquid  metal  to  the  pipe 
A  nai>  =  evaporator  area  =  ndopLe 


The  heat  transfer  to  the  Stirling  heater  is  used  by  the  Stirling  converter  to  produce 


electricity.  The  model  of  the  Stirling  converter  is  described  in  detail  in  the  next  section. 
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4.  Stirling  Converter  Model 

4.1  Stirling  Converter  Analysis 

Two  simulation  traits  that  are  involved  in  a  constant  tradeoff  are  simulation  speed  and 
sophistication.  While  a  main  goal  of  a  simulation  is  to  recreate  the  real  process  as 
accurately  as  possible,  accounting  for  all  the  variables  will  often  lead  to  simulations  that 
require  an  inordinate  amount  of  computer  resources  or  time  to  run.  The  Stirling  engine 
has  many  variables  and  simultaneous  processes  that  can  become  very  complex. 
Therefore,  in  order  to  develop  simulations  that  run  at  high  enough  speed  for  practical 
study  of  the  engine’s  performance  under  varying  conditions,  simplifying  assumptions 
were  made  to  reduce  the  level  of  complexity  of  the  model.  An  isothermal  analysis  is  one 
of  the  simplest  methods  for  developing  a  Stirling  engine  model  that  still  accounts  for  all 
of  the  system’s  dynamics. 

In  this  type  of  isothermal  analysis  it  is  assumed  that  temperatures  in  each  space  of  the 
engine  remain  constant.  By  making  this  assumption,  one  is  assuming  the  ideal 
performance  of  heat  exchangers  where  the  temperature  of  the  gas  is  assumed  equal  to  the 
heat  exchanger  wall.  The  efficacy  of  these  simplifications  is  evaluated  later  in  this 
section,  but  first  the  isothermal  method  is  explained  in  detail. 

The  isothermal  model  is  initially  derived  from  the  conservation  of  mass  in  the  working 
spaces  by  breaking  these  spaces  into  five  control  volumes:  the  compression,  cooler, 
regenerator,  heater,  and  expansion  spaces  (Figure  4-1).39’40 

mw  =  in  +  mk  +  mr  +  mh  +  me  4- 1 

where:  mw  =  total  mass  of  working  gas 

mc  =  compression  space  mass 
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nik  =  cooler  space  mass 
mr  =  regenerator  space  mass 
mi,  =  heater  space  mass 
me  =  expansion  space  mass 


Next,  the  assumption  is  made  that  the  pressure,  p ,  is  nearly  the  same  value  throughout 
the  working  spaces.  This  assumption  reflects  the  fact  that  the  difference  in  pressures 
between  the  working  spaces  is  small  and  can  be  accounted  for  later.  Using  the  isobaric 
assumption  and  the  ideal  gas  law,  the  following  equation  is  produced.39' 40  (Note  that  the 
same  subscripts  are  used  as  in  Equation  4-1) 


m  -  PK  - 

i_  pvk 

i_  pvr 

pvh 

|_  PVe 

w  R  T 

R  T, 

R  T 

R  T, 

R  T 

gas  c 

gas  k 

gas  r 

gas  h 

gas  e 

Since  Tc=Tk,  Equation  4-2  can  be  rewritten. 

m  ~  PV'  - 

i  Pyk 

a  pK 

Pyh 

,  pv. 

W  R  Ti 

gas  k 

R  T, 

gas  k 

R  T 

gas  r 

R  T, 

gas  h 

R  T, 

gas  h 

4-2 


4-3 


where:  p  =  working  space  pressure 

V  =  volume  of  the  designated  space 
Rgas  =  gas  constant  (helium) 

Ti  =  temperature  of  the  designated  space,  i 


Solving  for  the  working  space  pressure,  one  obtains  Equation  4-4. 
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P  = 


m»Rgns 


mwRgas 


Zl+Yl  +  Y^Yl  +  Y. 


4-4 


where: 


Y^+Yl+Y^+Yil+Y^ 

Tk  Tk  T  Th  Th 


Note  that  the  assumptions  made  allow  pressure  to  be  expressed  as  a  function  of  two 

on 

time  variant  parameters,  the  compression  and  expansion  space  volumes  (Vc  and  Ve). 
The  temperature  of  the  regenerator  is  assumed  to  be  the  effective  temperature  of  the 
regenerator  space  obtained  by  integrating  the  temperature  across  the  length  of  the 
regenerator.  Since  the  regenerator  wall  temperature  profile  is  produced  by  conduction 
from  the  heater  to  the  cooler,  the  temperature  of  the  wall  is  assumed  to  vary  linearly 


across  the  regenerator  length. 
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m. 


VrP  fP 

D  Jo 


dx 


VrP  \n(Th/Tk)  Vrp 


R 


gas 


KTh  -  Tk  )x  +  TkLr  ]  Rgas  (Th  -  Tk )  RgasTr 


T  = 


Wh/Tk) 

(Th~Tk) 


4-5 


4-6 


where:  L,  =  regenerator  length 

Tr(x)=(Th-Tk)x  +  TkLr 

=  temperature  profile  of  the  regenerator 
T,  =  effective  regenerator  temperature 

Next,  the  change  in  pressure  with  respect  to  time  is  calculated  by  taking  the  derivative 
of  Equation  4-5. 
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dp 


~mwRgas 


dV  dV. 


+ 


T  T 

1k  1h 


4-7 


The  rates  of  change  for  the  compression  and  expansion  space  volumes  found  in  these 
equations  were  then  determined  using  second  order  differential  equations  that  represent 
the  motion  of  the  piston  and  displacer.  These  differential  equations  can  be  put  into  the 
general  format  of  a  one  dimension  translational  mechanical  system.41 


mx  +  Dx  +  kx  =  F(t)  4-8 

where:  m  =  mass  of  displacer  or  piston 

I)  =  damping  coefficient 
ks  =  spring  constant 
x  =  mass  position 
x  =  mass  velocity 
x  =  mass  acceleration 
F(t)  =  forcing  function 

For  the  engine  configuration  used  in  this  research,  the  displacer  damping  force  is 
negligible  and  the  forcing  function  was  defined  by  the  pressures  in  the  bounce, 
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compression,  and  expansion  spaces  according  to  the  geometry  shown  in  Figure  4-1. 

The  displacement  value,  x,  of  zero  is  defined  for  the  equilibrium  position  and  the  positive 
direction  is  toward  the  expansion  space. 

Thus,  for  the  displacer  in  Figure  4-1,  the  damping  is  insignificant  and  the  equation  of 
motion  becomes, 


mdxd  +  kdxd  =  PhArod  +  PcApist  -  PeAd  4-9 

where:  Pb  =  bounce  space  pressure 

P(  =  compression  space  pressure 


*  In  order  to  keep  the  same  convention  as  previous  literature,  the  prefix  “d”  will  be  used  as  a  method  of 
labeling  the  derivative  with  respect  to  time. 
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Pe  =  expansion  space  pressure 
Arod  =  displacer  rod  cross  sectional  area 
APiSt  =  power  piston  cross  sectional  area 
Ad  =  displacer  cross  sectional  area 

The  power  piston  configuration  for  this  study  has  insignificant  spring  contributions 
and  its  damping  coefficient,  Dp,  is  mostly  a  function  of  the  alternator  and  is  defined  in 
detail  later  in  this  section.  The  equation  of  motion  for  the  power  piston  is, 

mpxp  +  Dpxp  =  (Pb  -  Pc)Ap  4-10 

Now  that  the  pressure  and  volume  variations  of  the  engine  are  known,  the  power  of  the 
engine  can  be  determined. 


clW  =  p(dVc  +dVe)=  p(Ap (xd  -xp)-Adxd)  4-11 

where:  dW  =  mechanical  power  from  the  working  gas 

dVc  =  time  derivative  of  compression  space  volume 
=  Ap(xd-xp) 

dVe  =  time  derivative  of  expansion  space  volume 

=  —Adxd 

Finally,  the  heat  transfer  and  energy  losses  are  all  that  need  to  be  determined  in  order 
to  complete  the  energy  balance  for  the  entire  engine.  Before  this  can  be  done,  however, 
the  mass  flow  rates  must  be  determined  for  each  of  the  spaces.  For  the  variable  volume 
compression  and  expansion  spaces,  this  is  done  using  the  differential  form  of  the  ideal 
gas  law  with  the  differential  temperature  term  defined  as  zero  (isothermal  assumption). 


dp  dV  dm  dT 
p  V  m  T 


dp  dV 
P  V 


4-12 


dm  =  m 


4-13 
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This  last  equation  can  be  simplified  by  creating  a  common  denominator  and  using  the 
ideal  gas  law. 


.  Vdp  +  pdV 

dm  =  m  - 

pv 

4-14 

Vdp  +  pdVc 
din  =  — - - 

V, 

4-15 

dm  =  V.dP+PdV. 

6  D  J1 

IX gas 1  h 

4-16 

The  heat  exchanger  mass  flow  rates  are  also  determined  using  the  differential  ideal  gas 
equation,  but  assuming  constant  volume. 


dp  dm 

4-17 

p  m 

dm  =  in— = 

Vdp 

4-18 

P 

R  T 

gas 

i  Vfdp 

dm,  =  — - - 

4-19 

R  T, 

gas  k 

dm,  =  V‘dp 

4-20 

h  R  T 

gas  n 

In  order  to  ensure  the  continuity  of  mass,  the  mass  flow  rates  at  the  boundaries 
between  the  spaces  must  be  defined.  A  positive  flow  is  defined  as  flow  towards  the 
expansion  space  in  all  of  the  following  relationships. 


drnck  =  —dmc  4-21 

dmkr  =  dmck  —  dmk  4-22 

dmhe  =  dme  4-23 

dmrh  =  dmhe  +  dmh  4-24 


where:  dmck  =  compression-cooler  boundary  mass  flow  rate 

dmkr  =  cooler-regenerator  boundary  mass  flow  rate 
dmhe  =  heater-expansion  boundary  mass  flow  rate 


dm,h  =  regenerator-heater  boundary  mass  flow  rate 
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The  boundary  mass  flows  for  the  regenerator  are  then  used  in  the  one-dimensional 
continuity  equation  to  determine  the  regenerator  mass  flow  rate. 

dmr  =  dmkr  —  dmrh  4-25 

Now  that  the  mass  flow  rates  are  known,  they  can  be  used  in  the  energy  equation  to 
determine  the  heat  transfer  of  the  heat  exchangers.  Ignoring  the  change  in  kinetic  energy 
and  the  energy  lost  to  friction,  the  energy  balance  for  a  generic  volume  can  be  calculated 
as  follows.40 


dQ  +  cp  (dmT  —  dmoTo )  =  dW  +  cvd  (ml')  4-26 

where:  dQ  =  heat  transfer  rate 

cp  =  average  specific  heat,  constant  pressure 
cv  =  average  specific  heat,  constant  volume 
dm,  =  mass  flow  rate  in 
dm0  -  mass  flow  rate  out 
dW  =  rate  of  work,  power 

For  isothermal  conditions,  this  equation  can  be  simplified  further.  The  constant 
volume  and  temperatures  of  the  heat  exchangers  results  in  the  appearance  of  no  heat 
transfer  from  these  spaces.  However,  the  heat  transfer  from  these  elements  is  what 
provides  the  energy  for  the  maintaining  the  temperature  of  the  variable  volume  spaces 
and  the  energy  needed  to  do  work.  Therefore,  we  can  redefine  our  control  volumes  as  the 
combination  of  the  heater  and  expansion  spaces  and  the  cooler  and  compression  spaces. 
Using  this  procedure,  the  following  relationship  is  determined.40 


f  dQ  +  RT  <f  dm  =  <f  dW 


4-27 


The  <fi  symbol  represents  the  integral  over  a  cycle.  For  steady  state,  the  net  mass  flow 

rate  is  zero  and  the  compression  and  expansion  work  rate  becomes  equal  to  the 
compression  and  expansion  heat  rate.40  The  subscript  i  is  used  to  denote  ideal  conditions. 
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dQk,  =  dWc  4-28 

dQhl  =  dWe  4-29 


From  this  relationship,  it  can  be  shown  that  the  efficiency  of  the  engine  is  the  Carnot 
efficiency.  Flowever,  to  determine  the  instantaneous  heat  transfer  and  to  determine  the 
heat  transfer  in  non-steady  state,  we  can  substitute  the  mass  flow  rate  into  Equation  4-25 


and  establish  the  following  relationship. 

a. 

P 

II 

1 

A, 

A 

4-30 

dQhi  =  -Vedp 

4-31 

In  order  to  overcome  the  artificiality  of  the  ideal  heat  exchangers,  the  hot  and  cold  end 
temperatures  were  assumed  to  be  the  value  necessary  to  provide  the  heat  transfer 
calculated  in  Equations  4-32  and  4-33.  This  will  decrease  the  temperature  range  that  the 


Cooler  Heater 


Figure  4-2.  Assumed  Wall  and  Gas  Temperature  Profiles  of  the  Stirling  converter. 
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engine  gas  experiences,  which  accounts  for  the  fact  that  the  heat  exchangers  are  not  ideal 
(Figure  4-2). 43 


T 

—  T  - 

dQhi 

4-32 

1h 

~  Jwh 

W vh 

T 

—  T  - 

dQki 

4-33 

1k 

~  Jwk 

Wvk 

where:  7),  =  heater  gas  temperature 

Tk  =  cooler  gas  temperature 
Twh  =  temperature  of  the  heater  wall 
Twk  =  temperature  of  the  cooler  wall 
dQhi  =  ideal  heater  heat  transfer 
dQki  =  ideal  cooler  heat  transfer 
hi,  =  heater  heat  transfer  coefficient 
h/(  =  cooler  heat  transfer  coefficient 
Awh  =  wetted  heater  area 
Awk  =  wetted  cooler  area 

The  regenerator’s  effectiveness  was  also  considered.  When  the  flow  goes  from  the 
heater  to  the  cooler,  the  flow  will  be  at  a  higher  temperature  than  the  cooler.  Also,  when 
the  flow  is  traveling  from  the  cooler  to  the  heater,  the  flow  temperature  will  be  lower  than 
that  of  the  heater.  In  both  of  these  cases,  the  heat  exchangers  will  have  a  higher 
magnitude  of  heat  transfer  than  the  original  value  calculated.  This  is  accounted  for  by 
calculating  the  regenerator’s  effectiveness.43 


dQk  =  dQki~\dQri\(l-e) 

4-34 

dQh  =  d Qhi  +  \dQrj  (1  -  e) 

4-35 

NTU 

£r  = 

'  1  i  A  JrTT  J 

4-36 

l  +  NTU 


where:  er  =  regenerator  effectiveness 

NTUr  =  number  of  transfer  units  for  the  regenerator 


h  A 

NTU  =-rfhsz. 
r  2Ar 


4-37 


52 


where:  h,  =  regenerator  heat  transfer  coefficient 

Awr  =  wetted  regenerator  area 
Ar  =  regenerator  cross  sectional  area 

The  two  in  the  denominator  of  Equation  4-37  is  due  to  the  fact  that  the  NTU  value  is 
accounting  for  both  the  heat  transfer  from  the  hot  gas  flow  to  the  regenerator  and  from 
the  regenerator  to  the  subsequent  cold  gas  flow.  Using  the  values  provided  for  the 
regenerator  from  NASA,  the  model  never  showed  effectiveness  values  below  0.98. 
Therefore,  disregarding  this  effect  would  not  introduce  significant  error  to  the  simulation. 

Due  to  the  isothermal  assumption,  the  heat  transfer  values  of  the  heat  exchangers  were 
taken  to  be  cycle  average  values.  To  determine  the  heat  transfer  coefficients,  the  mass 
flow  rates  of  the  spaces  were  used  to  find  the  flow  velocities.  From  the  definition  of 
mass  flow  rate,  the  velocity  of  the  gas  in  each  space  is  determined  as  follows. 


dmn  =  dmnV„ 

PnK  mnK 


4-38 


where:  un  =  flow  velocity  in  space  n 

pn  =  flow  density  in  space  n 

A„  =  cross-sectional  area  of  space  n 
n  =  subscript  designating  space 

The  mass  of  gas  in  each  space  is  defined  by  the  equation  of  state. 
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l  n  n 

R  T 

gas±  n 


4-39 


The  performance  of  this  engine  involves  mostly  laminar  flow  and  some  transitional 
flow.  For  the  heater  and  cooler,  which  both  consist  of  tubular  channels,  the  following 
correlations  were  used  for  calculating  the  friction  factor  and  Nusselt  number.44 
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for  Re  <  1502,/r  =  —  4-40 

Re 

for  Re>  1502,/r  =  4.41 

Re 

where:  Re  =  Reynolds  number 

=  pudh  / p 

p  =  fluid  density 

u  =  fluid  velocity 
dh  =  4/4/P,,: 

=  hydraulic  diameter 
Pw  =  wetted  perimeter 
p  =  fluid  viscosity 
fr  =  friction  factor 

for  Pe  <  1.5,  Nu  =  4.1 87(1 -0.0439Pe)  4-42 

for  Pe  >  1.5,  Nu  =  3.6568(1  +  1.227/  Pe2)  4-43 

where:  Pe  =  Peclet  number 

=  RePr 

Pr  =  Prandtl  number 

=  pcp/k 

k  =  fluid  thermal  conductivity 
Nu  =  Nusselt  number 

For  the  regenerator,  which  is  a  random  fiber  matrix  (Figure  4-3), 45  the  friction  factor, 
wetted  perimeter,  and  Nusselt  number  are  calculated  as  follows.46 


fr  =  —  +  4.53Re~a67 
Re 
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Fig  4-3.  Random  fiber  regenerator  matrices  with  80%,  left,  and  88%,  center,  porosity  and 
their  location  within  a  Stirling  converter,  (site) 
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4(1  -f3)  A 

W  7 

dw 

Nu  =  l  +  l.l6Pe0Mp2M 


4-45 

4-46 


where:  [>  =  regenerator  porosity 

=  void  volume/total  volume 
dw  =  wire  diameter 

The  heat  transfer  coefficients  for  the  Stirling  heat  exchangers  can  then  be  calculated 
using  the  Nusselt  numbers,  fluid  thermal  conductivity,  and  the  hydraulic  diameter  of  the 
space. 


h 


Nu-k 

dh 
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In  order  to  determine  the  pressure  in  the  expansion  space,  the  flow  losses  in  the 
Stirling  heat  exchangers  are  determined  using  the  friction  factors  calculated  in  Equations 
4-39,  4-40,  and  4-43.  These  flow  losses  correlate  to  a  change  in  pressure  across  each  of 
the  heat  exchangers  (cooler,  heater,  and  regenerator)  in  the  Stirling  engine. 
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where:  A p  =  pressure  difference  across  a  heat  exchanger 

F  =  frictional  force  in  a  heat  exchanger 
Ahx  =  cross  sectional  flow  area  of  a  heat  exchanger 


F  =  -v\fL+^\^k 

v  dh  Lhx  2 
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where:  V  =  heat  exchanger  fluid  volume 

K  =  local  loss  coefficient 
Lhx  =  heat  exchanger  length 


The  local  loss  coefficient,  K,  accounts  for  bends  in  the  flow  and  is  usually  assumed  to 
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be  1.5  for  heater  and  cooler  due  to  the  ninety  degree  bends  into  the  compression  and 
expansion  spaces  (Figure  4-1). 

P  =  p  4-50 

Pe  =  Pc  +  APk  +  APr  +  APh  4-5  1 

where:  Pc  =  compression  space  pressure 

Pe  =  expansion  space  pressure 
A/p  =  cooler  pressure  loss 
A pr  =  regenerator  pressure  loss 
A pi,  =  heater  pressure  loss 

4.2  Alternator  Analysis 

With  the  pressures  known,  the  last  parameter  that  must  be  calculated  to  determine  the 
engine  dynamics  is  the  piston  damping  due  to  the  alternator.  The  piston-alternator 
configuration  consists  of  the  piston,  magnets  that  are  attached  to  the  piston,  and  wire  coils 
that  suiTound  the  piston  magnets.  The  magnets  are  very  close  to  the  coils.  This  allows 
the  alternator  voltage  to  be  calculated  by  a  linear  relationship  to  the  piston  velocity  as 
follows: 


E  =  Kit*,  4-52 

where:  E  =  alternator  voltage 

kait  =  alternator  constant  (38.5  Ys/m) 
x  =  piston  velocity 

This  voltage  induces  a  current  in  the  coils  that  must  be  calculated  using  a  second  order 
equation  describing  the  electrical  circuit.  The  model  currently  does  not  have  a  capacitor, 
but  capacitance  is  included  in  the  following  equation  since  tuning  capacitors  are  often 


found  in  Stirling  converters. 
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E  =  *„„ji„=L/  +  «  +  U  4-53 

where:  E  =  time  derivative  of  the  voltage 

L  =  alternator  inductance 
R  =  total  resistive  load 

=  alternator  +  load  resistance 
C  =  circuit  capacitance  (ignored) 

I  =  second  time  derivative  of  the  current 
1  =  first  time  derivative  of  the  current 
/  =  circuit  current 


Figure  4-4.  Schematic  of  the  alternator  circuit 
Using  the  current  calculated  in  Equation  4-53,  the  electrical  power  produced  by  an 
engine  was  calculated  by  multiplying  the  current,  voltage,  and  efficiency. 


We=rialtIE 
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where:  We  =  electrical  power  produced  by  a  single  engine  (50  kW(e)) 

rjait  =  alternator  efficiency 

Since  there  are  four  engines  in  the  system,  the  value  calculated  in  Equation  4-54  must 
be  multiplied  by  four  to  determine  the  total  electrical  energy  produced.  The  net 
efficiency  of  the  Stirling  converter  is  then  calculated  by  comparing  the  electrical  energy 
produced  to  the  thermal  energy  added. 


V 


cony 
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where: 


rjconv  =  instantaneous  converter  efficiency 


57 


The  average  converter  efficiency  as  a  function  of  time  can  be  determined  by 
integrating  the  electrical  power  over  time  and  dividing  the  result  by  the  integral  of  the 
thermal  power  added  to  the  converter. 

The  damping  force  on  the  piston  was  also  calculated  by  multiplying  the  current  value 
from  Equation  4-53  by  the  alternator  constant. 


Fp=- 


-KnI : 


-k~  x 

alt  P 


-D  x 

p  p 


■‘alt 


D 


p 
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where:  Fp  =  damping  force  on  the  piston 

Zait  =  total  impedance  of  the  alternator  circuit 
Dp  =  damping  coefficient 

The  total  impedance  includes  the  resistive  load  and  the  impedance  of  the  inductor. 


z.„=l  R  +  Zj  =  J¥+cFfLf  4-58 


where:  Zinci  =  inductor  impedance 

=  2?r/Lj 
/=  frequency 

The  total  electrical  impedance  was  determined  by  assuming  sinusoidal  displacement 
and  velocity  of  the  piston  (Equations  4-59  and  4-62). 

xp  =  Xp  sin (2-71/0  4-59 


where  xp  =  piston  displacement 

Xp  =  maximum  piston  displacement  (13  mm) 
/=  frequency  (96.1  Hz) 
t  =  time 


The  frequency  of  the  piston  was  determined  by  using  the  second  order  system  Laplace 

42 


transfer  function. 
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X(s) 


C 0 


4-60 


F(s )  s  +  2C,us  +  or 

where:  co  =  angular  frequency 

=  2nf 

C  =  damping  ratio 

Using  Equation  4-8,  it  can  be  shown  the  following  relationships  exist. 


/=_!_  L 
2tt  V  m 
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To  determine  the  frequency,  the  values  for  the  displacer  spring  constant  (kd=3.1xl06) 
and  mass  (m=8.5  kg)  were  used  in  Equation  4-59 


_ L  3Ae 6 
/_27rV  8.5 
=  96.1  Hz 

By  taking  the  derivative  of  Equation  4-59,  the  piston’s  velocity  was  approximated  as 
follows. 


xp  =  Zp(2tt/)cos(2tt/0  4-62 


where:  x  =  piston  velocity 

Then,  Equations  4-52  and  4-62  were  substituted  into  the  electrical  power  relationship 
(Equation  4-63)  and  this  equation  was  solved  for  the  total  electrical  impedance  of  the 
alternator  circuit. 


E 

W=-  p 


2  Z 


alt 


{kahxpYp  (2m f)zK,tX; 

2ZaIt 


2  Z 


alt 


Ep  =  peak  voltage 
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Equation  4-63  was  then  solved  for  the  total  impedance,  resulting  in  Equation  4-64. 
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(2-n ffk2altXl 

2We 

(2tt-96.1)2(38.5)2(0.013)2 

2(50000) 


Z„/,  =  0.914  O 
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The  load  resistance  was  determined  by  using  a  steady-state  relationship  to 
represent  the  impedance  of  the  inductor. 

Zah  =  \J(Rload  +Ralt)2  +  Zld  4'65 

where:  Rioad  =  load  resistance 

Rah  =  alternator  resistance  (0.028  Q) 

Zj„d  =  alternator  inductor  impedance 

Rload  =  \jZalt  ~~  Zind  ~~  Ralt 

where:  Lind  =  alternator  inductance  (0.2  mH) 

Rload  =  >/(0.9 14)2  -  (2ir96.  l(2e  -  4))2  -  0.028 
=  0.878  n 


Initial  tests  of  the  system  produced  errors.  These  errors  were  corrected  by  modeling 
the  motion  of  the  piston  and  displacer  as  perfect  sinusoids  with  the  phase  angle  given  by 
NASA.  This  method  was  used  instead  of  integrating  the  acceleration  determined  from 
Newton’s  Second  Law.  This  sinusoidal  simulation  produced  negative  volumes,  so  the 
volume  of  the  compression  space  was  then  increased  until  all  volumes  remained  positive 
and  the  minimum  volumes  were  no  less  than  five  percent  of  their  respective  mean  values. 
This  was  achieved  when  the  compression  space  value  was  increased  thirty  percent  above 
the  initial  value  provided  by  NASA. 


Next,  the  sinusoidal  restraint  was  removed  from  the  piston  and  the  simulation  was  run 
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again  using  Newton’s  Second  Law.  During  this  simulation,  the  motion  of  the  piston 
showed  an  overdamped  response.  This  issue  was  resolved  by  evaluating  the  differential 
equations  of  motion  (force  balances)  related  to  the  power  piston. 


mpxp  +  K,,1  =  (pb  -p  )A»  =  F„ (0 


E  =  Ki,xP  =  LI +R1 

m  L  ..  m  R 

FJt)  = 


I  +  KJ 


xalt 


xalt 
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Similarly  to  Equation  4-59,  the  mass  of  the  power  piston  (mp)  can  be  related  to  its 


frequency  (f  or  w)  by  determining  the  transfer  function  of  equation  4-65. 


38. 52 

_  (27r96.1)2(2e-4) 

=  20.3  kg 

When  this  value  of  the  power  piston  mass  was  substituted  for  the  value  provided  by 
NASA,  25  kilograms,  the  Stirling  model  reached  steady  state  and  closely  matched  the 
expected  values  for  piston  motion,  displacer  motion,  and  electrical  power.  The  results  of 
this  simulation  are  shown  in  section  6  with  the  other  simulation  results. 

Using  all  of  these  correlations,  a  simulation  in  MATLAB/Simulink  was  created  to 
determine  the  performance  of  the  Stirling  engine.  The  analysis  of  these  simulations  is 
found  in  the  Section  6.  The  listing  of  the  engine  parameters  provided  is  attached  in 
Appendix  E.42 
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5.  Heat  Rejection  System 

The  liquid  metal  heat  rejection  loop  feeds  through  the  coolers  of  the  Stirling  converters 
to  absorb  the  waste  heat  and  transport  it  to  the  radiators  via  heat  pipes  using  NaK  as  a 
working  fluid.  NaK  has  large  liquid  temperature  range  (260-1060K  under  standard 
pressure),  which  allows  it  to  be  used  for  each  of  the  three  pumped  loops  of  this  system. 
The  arrangement  of  the  heat  rejection  system  is  displayed  in  Figure  5-1. 


NaK  Return  Flow 

Figure  5-1.  Heat  Rejection  System  Configuration  (Not  to  scale  -  system  has  36  panels  vice  10). 


Figure  5-2.  The  side  and  top  view  of  the  interface  of  the  cooler  and  the  heat  rejection  loop. 


The  geometry  of  the  interface  of  this  pumped  loop  with  a  cooler  is  shown  in  Figure  5- 
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2.  The  liquid  metal  flows  through  channels  that  run  concentrically  around  the  rows  of 
holes  that  the  gas  (helium)  in  the  cooler  flows  through.  A  figure  of  the  Stirling  converter 
cooler  is  given  in  Figure  4- 1 .  The  heat  transfer  to  the  liquid  metal  was  treated  as  flow 
between  two  parallel  plates,  since  this  correlation  was  the  most  similar  to  the  geometry  of 
the  channels. 


^ chan 


KNaK 

du 


Nu 


chan 


25 


0.00197 
=  1.203xl05 


-9.49 
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where:  Khan  =  channel  heat  transfer  coefficient 

kNaK  =  NaK  thermal  conductivity 
dh  =  channel  hydraulic  diameter 
Nuchan  =  channel  flow  Nusselt  number  (9.49) 


The  change  in  temperature  of  the  liquid  metal  in  the  Stirling  engine  cooler  channels 
was  calculated  using  the  energy  equation. 


dT 


chan 


1 


dt 


CpPNaK^c 


_(h  A  (T  —T  V 

l ' 1 chan 1  wchan  V ±  wk  ±  chan  ’ 


'  C pd1  chan  (Tarl 


chan 


where:  Tchan  =  channel  coolant  temperature 

cp  =  coolant  specific  heat  (NaK) 

PNaK  =  coolant  density  (NaK) 

Vchan  =  channel  volume 
A  wchan  =  channel  wetted  area 
Twk  =  cooler  wall  temperature 
mchan  =  coolant  mass  flow  rate 

Tari  =  channel  exit  coolant  temperature 
Txr,  =  channel  entrance  coolant  temperature 

From  the  cooler,  the  liquid  metal  flows  toward  the  radiators  through  a  circular  pipe. 

Heat  pipes  are  used  to  transfer  heat  from  the  pumped  loop  flow  to  the  radiator  panels. 


The  evaporator  end  of  the  heat  pipe  (see  Figure  2-5)  is  embedded  into  the  transfer  loop 
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pipe  so  that  it  is  perpendicular  to  the  oncoming  flow.  As  shown  in  Figure  5-1,  two  heat 
pipes  were  placed  at  each  location  in  order  to  create  symmetry  for  the  radiator  panels  on 
both  sides  of  the  coolant  pipe.  The  Nusselt  number  for  this  liquid  metal  flow  geometry 
was  calculated  using  Equation  5-2. 32 

^=5  +  0.0257*“  5-3 


where:  Nurhpe  =  Nusselt  number  of  radiator  heat  pipe 

evaporator 

Perhpe  =  flow  Peclet  number  at  the  heat  pipe 
The  condenser  end  of  the  heat  pipe  is  embedded  into  the  radiator  panel.  The  heat  pipes 
used  for  this  portion  of  the  system  consist  of  titanium  walls  and  wicks  with  water  as  the 
working  fluid.  The  heat  transfer  through  the  heat  pipes  was  modeled  with  thermal 
resistances  as  described  in  Section  3.4.  The  temperature  of  the  NaK  coolant  at  each 
section  of  two  heat  pipes  is  calculated  similarly  to  Equation  5-2. 


dt 
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oPNaK^cs 


T  —T 

ab  abp 


R , 


hptot 
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where:  Tab  =  section  average  coolant  temperature 

Vcs  =  section  coolant  volume 
mrl  =  heat  rejection  loop  mass  flow  rate 
Tari  =  section  inlet  temperature 
Tbri  =  section  outlet  temperature 
Tabp  =  section  radiator  panel  temperature 
R hptot  =  radiator  heat  pipe  total  thermal  resistance 

Note  that  the  section  mass  flow  rate  of  the  NaK  in  the  rejection  loop  is  four  times  larger 


than  the  channel  mass  flow  rate  since  there  are  four  Stirling  converters  in  the  model. 


Also,  the  section  inlet  temperature  is  assumed  equal  to  the  outlet  temperature  of  the 


previous  section,  except  after  the  last  radiator  where  a  time  delay  is  added  to  account  for 


the  time  necessary  for  the  flow  to  return  to  the  Stirling  converters.  Each  interval  between 


64 


sections  is  labeled  alphabetically  from  a  to  s  (Figure  5-1),  except  for  the  inlet  temperature 
at  the  Stirling  cooler. 

Since  titanium  is  a  poor  conductor,  a  graphite  foam  interface  was  used  to  connect  the 
heat  pipe  condenser  to  the  radiating  fins  of  the  radiator  panels  (Figure  5-3).  This 


XHeat  Pipe  Envelope^  ^\GFrC  Fins 

Figure  5-3.  Diagram  of  radiator  panel  geometry  [47]. 

configuration  comes  from  recent  research  completed  at  Glenn  Research  Center.47  The 
radiating  fins  were  made  of  graphite  fiber  reinforced  composite  (GFRC),  and  were 
assumed  to  have  an  emissivity  of  O.9.47  Since  the  foam  interface  and  radiator  fins  both 
primarily  consist  of  graphite  and  have  high  conductivity  perpendicular  to  the  heat  pipes, 
they  were  treated  as  a  single  thermal  mass  of  the  same  temperature.  The  effect  of  heat 
transfer  through  the  aluminum  honeycomb  is  neglected  (Figure  5-3).  Due  to  symmetry, 
the  fin  temperature  for  each  radiator  pair  was  assumed  to  be  the  same.  An  energy 
balance  was  used  to  calculate  the  temperature  of  each  pair  of  panels  as  given  in  Equation 


5-5. 
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where:  cpr  =  radiator  specific  heat 

pgr  =  graphite  density 

Vmd  =  radiator  volume 
e  =  radiator  emissivity 
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g  =  Stephan-Boltzman  constant  =5.67E-8  W/(m2K4) 

Amd  =  radiator  surface  area  =  4  x  radiator  fin  area 
Tamb  =  ambient  space  temperature  (usually  4  K) 

The  area  of  each  set  of  radiator  panels  was  based  the  width  provided  by  a  study  of 
identical  system  and  a  length  perpendicular  to  the  pumped  loop  such  that  the  panels 
remained  inside  a  10  degree  cone  that  is  protected  by  shielding.  The  final  radiator 
model  had  the  characteristics  given  in  Table  5-1. 


Variable 

Value 

Number  of  radiator  panels 

36 

Total  radiator  area  (m2) 

404.9 

Total  radiator  mass  (kg) 

1,313 

Radiator  heat  pipe  diameter  (cm) 

2.54 

Heat  pipe  wick  outer  diameter  (cm) 

2.29 

Heat  pipe  wick  inner  diameter  (cm) 

2.04 

Heat  rejection  loop  mass  flow  rate  (kg/s) 

13.7 

Rejection  loop  return  time  delay  (s) 

2.5 

Steady  state  heat  radiated  (kW) 

415 

Table  5-1.  Radiator  model  parameters. 


The  simulation  results  of  the  radiator  model  can  be  found  in  Section  6. 
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6.  Simulation  Results 

This  section  describes  the  simulation  results  obtained  from  Simulink.  Initially,  each 
of  the  three  main  system  components  was  modeled  individually  and  tested  (Sections  6.1 
through  6.3).  After  each  sub-model  was  verified,  a  model  of  the  entire  system  was 
created  by  integrating  each  of  the  sub-models.  Simulations  of  the  entire  system’s 
response  to  various  perturbations  were  then  run  (Section  6.4). 


6.1  Partial  Model  Simulation  Results  for  the  Reactor,  Heat  Exchanger, 
and  Transfer  Loops 


After  the  completion  of  the  individual  Simulink  models  of  the  reactor,  heat  transfer 
loops,  and  the  heat  exchanger  that  links  the  transfer  loops,  the  individual  models  were 
linked  and  tested  with  a  constant  heat  sink  temperature.  At  steady  state,  this  model 
produced  the  parameters  given  in  Table  6-1. 


Sub-System  Parameter 

Steady  State  Value 

Reactor  Power,  kW(t) 

600 

Reactor  Precursor  Power  Equivalence,  GW  (Gigawatt) 

50.6 

Reactor  Reactivity 

0 

Fuel  Temperature,  K 

1097.6 

Average  Reactor  Coolant  Temperature,  K 

950 

Reactor  Coolant  Inlet  Temperature,  K 

900 

Reactor  Coolant  Outlet  Temperature,  K 

1000 

Heat  Exchanger  Primary  Coolant  Inlet  Temperature,  K 

1000 

Heat  Exchanger  Primary  Coolant  Outlet  Temperature,  K 

900 

Heat  Exchanger  Primary  Coolant  Average  Temperature,  K 

950 

Heat  Exchanger  Primary  Wall  Temperature,  K 

943.7 

Heat  Exchanger  Secondary  Wall  Temperature,  K 

932.7 

Heat  Exchanger  Secondary  Coolant  Inlet  Temperature,  K 

870.7 

Heat  Exchanger  Secondary  Coolant  Outlet  Temperature,  K 

970.7 

Heat  Exchanger  Secondary  Coolant  Average  Temperature,  K 

920.7 

Average  Secondary  Coolant  Temperature  at  Sink,  K 

920.7 

Secondary  Heat  Transfer  Loop  Sink  Temperature,  K 

908.8 

Table  6-1.  Steady  State  Values  for  the  Reactor-Heat  Exchanger-Transfer  Loops  Partial  Model. 
The  results  supported  the  validity  of  this  portion  of  the  model  by  demonstrating  the 

100  K  temperature  difference  across  the  reactor  that  was  previously  produced  in  similar 


simulations  by  Naval  Reactors.  The  temperature  drop  across  the  heat  exchanger  and 
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secondary  loop  heat  sink  were  also  100  K  during  the  simulation.  This  shows  that  the  heat 
transfer  values  across  the  system  were  all  the  same,  as  expected  for  steady  state 
conditions. 

For  the  first  dynamic  simulation  (Table  6-2),  the  heat  sink  (Stirling  heater)  for  the 
secondary  transfer  loop  was  given  a  constant  heat  rejection  of  600  kW(t)  while  the 
reactivity  of  the  reactor  was  given  a  step  input  of  0.001  at  5  seconds.  This  simulation  is 
meant  to  model  the  effects  of  pulling  a  control  rod  partially  out  of  the  reactor.  The 
reactor  response  for  this  transient  is  shown  in  Figures  6-1  and  6-2  (page  68).  The 
responses  of  the  heat  exchanger  and  secondary  loop  are  shown  in  Figures  6-3  and  6-4 
(page  69),  respectively.  Further  model  results  for  this  case  are  provided  in  Appendix  A. 


Reactor  Test 

Simulation 

Variable  Changed 

1 

0.001  Reactivity  Increase 

Reactivity 

2 

-0.001  Reactivity  Decrease 

Reactivity 

3 

10  K  Increase  in  Stirling  Heater  Temperature 

Heater  Temperature 

4 

10  K  Decrease  in  Stirling  Heater  Temperature 

Heater  Temperature 

Table  6-2.  Dynamic  Simulations  for  Partial  Model  of  Reactor,  Heat  Exchanger,  and  Transfer  Loops. 

These  initial  simulation  results  support  that  the  model  is  working  properly.  As 

expected,  the  positive  reactivity  insertion  increased  the  reactor  power  level.  This  increase 
in  power  level  caused  an  immediate,  sharp  increase  in  the  fuel  temperature  (~  40  K  in  one 
minute).  The  NaK  coolant  temperatures  increase  shortly  thereafter.  The  outlet 
temperature  rose  the  fastest,  while  the  inlet  temperature  had  a  delayed  response  due  to  the 
transfer  loop.  The  temperature  increase  propagated  throughout  the  rest  of  the  transfer 
loops  as  expected. 

Due  to  the  negative  temperature  coefficient  of  the  reactor,  the  increasing  temperature 
in  the  reactor  after  the  positive  reactivity  insertion  caused  natural  feedback  effects,  which 


returned  the  reactivity  back  to  zero  and  the  reactor  to  its  original  power  level.  Also 


x  10" 


Figure  6-1.  Partial  model  reactor  response  to  0.001  reactivity  step  input  at  5  seconds.  The  reactor  power 
increased  with  the  reactivity  increase  and  returned  to  the  original  value  after  395  seconds. 


Figure  6-2.  Partial  model  reactor  response  to  0.001  reactivity  step  input  at  5  seconds.  The  reactor 
temperatures  increase  as  a  result,  with  fuel  reacting  the  fastest  and  the  inlet  temperature  reacting  the 

slowest. 

shown  in  Figure  6-1  is  the  heat  flow  from  the  fuel  to  the  coolant  and  the  precursor  power 
equivalence,  which  takes  into  account  the  affects  of  the  fission  products. 

When  analyzed  as  a  whole,  the  sub-system  in  the  partial  model  had  a  constant  heat 


rejection  rate,  but  the  heat  addition  increased  with  the  reactivity  insertion  and  returned 


Figure  6-3.  Partial  model  heat  exchanger  response  to  0.001  reactivity  step  input  at  5 
seconds.  The  temperatures  of  the  heat  exchanger  increase  with  a  delay. 


Figure  6-4.  Partial  model  secondary  loop  response  to  0.001  reactivity  step  input  at  5 
seconds.  The  temperatures  of  the  secondary  loop  at  the  Stirling  Heater  increase  with  an 

back  to  the  original  value  after  about  395  seconds.  Therefore,  the  system  had  a  net  heat 
addition  throughout  the  simulation,  which  explains  why  the  final  temperatures  were 
higher  than  the  initial  temperatures.  This  test  showed  that  this  partial  model  functions 


correctly  and  that  these  components  would  tend  to  be  stable  and  load  following. 
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A  similar  simulation  was  run  with  a  negative  step  input  of  -0.001  reactivity  (Appendix 
B).  This  test  simulates  the  partial  insertion  of  a  control  rod.  The  system  responded  with 
the  same  magnitudes  as  the  0.001  reactivity  insertion,  but  in  the  opposite  direction.  This 
test  also  showed  that  this  partial  model  functions  correctly  and  that  these  components 
would  tend  to  be  stable  and  load  following. 

The  third  dynamic  simulation  involved  decreasing  the  secondary  heat  transfer  loop 
sink  temperature  with  a  negative  one  degree  Kelvin  per  second  ramp  input  for  ten 
seconds  and  observing  the  system  response  (Appendix  C).  This  simulation  replicates  a 
increased  power  draw  from  the  Stirling  converters.  The  response  for  this  case  is  shown 
in  Figures  6-5  through  6-8  (pages  71-72).  Figure  6-5  shows  that  the  reactor  power 
increased  to  meet  the  increased  heat  removal  that  occurred  as  the  heat  sink  temperature 
decreased  and  that  this  response  was  delayed  due  to  the  transit  time  of  the  NaK  in  the 
transfer  loops.  However,  this  increase  in  power  was  small  (only  1%  of  the  original  value) 
due  to  the  small  magnitude  of  the  perturbation  and  the  relatively  fast  response  of  the 
system  (the  reactor  begins  its  response  5  seconds  after  the  perturbation).  The  sub-system 
overcompensates  slightly  since  the  heat  rejected  by  the  second  loop  equalizes  (after  ~  195 
seconds)  with  the  heat  addition  from  the  reactor  at  a  value  slightly  higher  than  the 
original  value.  This  is  verified  by  the  final  reactor  power,  reactivity,  and  heat  sink 
temperature  values.  The  results  of  test  three  provide  further  support  that  the  sub-system 
in  this  partial  model  is  stable  and  tends  to  be  load  following.  The  extra  heat  drawn  from 
the  Stirling  converters  was  met  with  an  increase  in  power  in  the  reactor. 

Test  four  conducted  the  same  test  but  with  a  positive  temperature  increase  (Appendix 
D).  The  simulation  had  the  same  response  times  and  magnitudes,  but  in  the  opposite 
direction.  Test  four  provides  further  evidence  that  the  system  is  stable  and  load 
following. 
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Figure  6-5.  Reactor  response  to  decreasing  the  secondary  transfer  loop  sink  temperature  1 
Kelvin/second  for  10  seconds.  The  reactivity  and  power  of  the  reactor  increase  as  colder  NaK  enters  the 
reactor  and  they  decrease  as  the  NaK  is  heated  above  the  equilibrium  value. 
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Figure  6-6.  Reactor  response  to  decreasing  the  secondary  transfer  loop  sink  temperature  1 
Kelvin/second  for  10  seconds.  The  coolant  temperatures  drop,  resulting  in  a  subsequent  temperature 

drop  in  the  fuel. 
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Figure  6-7.  Primary  side  heat  exchanger  response  to  decreasing  the  secondary  transfer  loop  sink 
temperature  1  Kelvin/second  for  10  seconds.  The  outlet  temperature  responds  earlier  than  the  inlet 
temperature  due  to  the  transit  time  of  NaK  in  the  primary  loop. 
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Figure  6-8.  Secondary  side  heat  exchanger  response  to  decreasing  the  secondary  transfer  loop  sink 
temperature  1  Kelvin/second  for  10  seconds.  The  decrease  in  the  inlet  temperature  draws  more  heat  from 
the  primary  loop,  resulting  in  the  increase  in  outlet  temperature.  Both  temperatures  return  to  normal 
after  the  primary  coolant  temperature  returns  to  its  original  value. 


73 


6.2  Partial  Model  Simulation  Results  for  a  Stirling  Converter 

A  Stirling  model  was  created  using  the  method  outlined  in  Section  4.  Data  provided 
by  Glenn  Research  Center  (GRC)  was  used  as  a  guide  (Appendix  E).  The  model  was 
tested  by  analyzing  the  steady  state  response  when  a  constant  heat  input  was  added  to  the 
Stirling  heater.  The  initial  0.3  seconds  of  this  response  is  shown  in  Figures  6-9  and  6-10 
(page  74).  It  is  clear  from  counting  the  oscillations  in  these  graphs  that  the  Stirling  is 
operating  near  the  expected  frequency  of  96  Hz  (slightly  higher  than  the  90  Hz  provided 
by  GRC).  As  demonstrated  in  the  left  plot  of  Figure  6-9,  the  expansion  and  compression 
space  temperatures  are  nearly  identical.  This  validates  the  small  pressure  difference 
assumption  made  in  Section  4-1.  The  magnitude  of  the  pistons  movement  is  just  short  of 
the  expected  13  mm  and  the  mean  pressure  is  near  the  expected  13.6  MPa.  Also,  taking 
the  area  of  the  pressure- volume  diagram  (Figure  6-10)  and  multiplying  by  the  frequency 
and  electrical  efficiency  produces  45  kW(e).  This  value  matches  the  electrical  power 
provided  by  the  simulation  (Figure  6-11,  page  75).  This  value  can  be  obtained  by 
subtracting  the  heat  rejected  from  the  heat  added  in  Figure  6-11  and  multiplying  by  the 
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Figure  6-9.  Stirling  steady  state  response.  The  values  initially  varied  as  the  model  was  started. 
This  start-up  transient  ended  within  0.13  seconds. 
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alternator  electrical  efficiency  (93%).  All  of  the  preceding  evidence  proves  that  the 
Stirling  model  is  working  correctly  and  the  simulation  reveals  that  the  model  is  stable.  It 
is  worth  noting,  however,  that  the  simulations  shown  in  this  section  are  for  one  Stirling 
converter,  while  the  completed  system  has  four  converters.  Also,  the  Stirling  electrical 
output  is  lower  than  the  50  kW(e)  needed  to  produce  200  kW(e).  This  issue  was  resolved 
when  the  complete  system  model  was  created  and  it  will  be  discussed  in  further  detail  in 
Section  6.4.  The  steady  state  values  determined  by  the  model  are  given  in  Table  6-3. 


Sub-System  Parameter 

Steady  State  Value 

Electrical  Power,  kW(e) 

45 

Heat  Addition,  kW(t) 

143 

Heat  Rejection,  kW(t) 

94 

Frequency,  Hz 

96 

Mean  Pressure,  MPa 

13.6 

Piston  Amplitude,  mm 

12.5 

Displacer  Amplitude,  mm 

6.5 

Heater  Temperature,  K 

925 

Cooler  Temperature,  K 

463 

Alternator  Constant,  Vs/m 

38.5 

Alternator  inductance,  mH 

0.2 

Alternator  resistance,  Q 

0.028 

Load  resistance,  Q 

0.878 

Table  6-3.  Steady  state  values  for  the  Stirling  converter  sub-model. 
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Figure  6-10.  Steady  state  plots  of  the  Stirling  PV  diagram  and  Volume  versus  time. 


Figure  6-11.  Stirling  power  values  for  a  steady  state  response.  The  values  initially  varied  as  the 
model  was  started.  This  start-up  transient  ended  within  0.13  seconds. 
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6.3  Heat  Rejection  System  Model  Simulation  Results 
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The  model  of  the  heat  rejection  system  described  in  Section  5  was  also  modeled  and 
tested  in  Simulink.  The  steady  state  temperatures  of  the  model  agreed  with  theory.  In 
order  to  simulate  conditions  similar  to  a  reactivity  insertion,  the  thermal  power  into  the 
NaK  channels  of  the  Stirling  coolers  were  increased  from  the  steady  state  value  by  20 


Flow  and  Radiator  Panel  Temperatures  of  Heat  Rejection  System  for  a 
20%  Increase  in  Thermal  Input  at  1000  s 


Figure  6-12.  Transient  Test  of  the  Radiator  Model 
percent  (600  to  720  kW(t)).  The  results  of  the  transient  test  are  shown  in  Figure  6-12.  It 

is  clear  from  this  figure  that  the  NaK  channel  flow  has  a  smaller  thermal  capacitance 

since  it  increases  in  temperature  so  quickly  (solid  lines).  Also,  one  can  see  the  increase  in 

thermal  capacitance  of  the  radiators  as  distance  from  the  heat  source  increases  along  the 

pumped  loop  and  the  radiator  size  increases  (the  17th  set  of  panels  heats  up  slower  than 

the  1st  set).  This  simulation  suggests  that  the  heat  rejection  model  is  working  correctly 


and  that  the  model  is  stable.  The  steady  state  values  calculated  by  this  model  are  given  in 
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Table  6-4. 

After  the  validation  of  the  radiator  model,  it  was  linked  with  the  Stirling  and  reactor 
models  to  create  the  model  of  the  entire  system.  The  simulation  tests  of  the  system 
model  are  displayed  and  analyzed  in  the  next  section. 


Sub-System  Parameter 

Steady  State  Value 

Number  of  radiator  panels 

36 

Number  of  heat  pipes 

36 

NaK  temperature  at  1st  set  of  heat  pipes,  K 

457 

NaK  temperature  at  8th  set  of  heat  pipes,  K 

449 

NaK  temperature  at  1 7th  set  of  heat  pipes,  K 

428 

Temperature  of  1st  set  of  radiator  panels,  K 

415 

Temperature  of  8th  set  of  radiator  panels,  K 

387 

Temperature  of  17th  set  of  radiator  panels,  K 

355 

Heat  Rejection,  kW(t) 

600 

Table  6-4.  Steady  state  values  of  the  Heat  Rejection  System. 
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6.4  System  Simulation  Results 

After  the  system  model  was  created,  tests  (Table  6-5)  were  run  to  determine  that  the 
component  models  were  linked  together  properly.  This  was  accomplished  by 
determining  the  initial  conditions  of  the  system  analytically  and  then  running  the  model 
and  comparing  the  results  to  the  expected  results.  Each  simulation,  including  all  of  the 
previous  simulations,  were  run  using  the  variable-step  ode45  (Dormand-Prince)  as  the 
numerical  solver  with  a  tolerance  of  1 0  4.  This  solver  is  a  Runga-Kutta  solver  and  is  the 
most  accurate  option  available  in  Simulink. 

Initially,  each  model  did  not  integrate  smoothly  due  to  small  differences  in  the  initial 
conditions  created  by  using  separate  reference  points.  This  problem  was  solved  by  taking 
the  values  of  the  Stirling  heater  and  cooler  temperatures  and  using  them  as  a  reference  to 
determine  the  initial  conditions  of  the  reactor  and  radiator  models.  As  mentioned  in 
Section  6.2,  the  four  Stirling  converters  were  not  producing  the  required  electrical  power 
(200  kW(e)),  because  their  efficiency  was  slightly  lower  than  that  initially  assumed 
(31.7%  vice  33.3%).  The  effect  of  this  lower  efficiency  was  that  reactor  thermal  power 
needed  to  be  slightly  greater  than  600  kW(e)  as  shown  in  Figure  6-13. 

As  can  be  seen,  the  steady  state  condition  of  this  system  nearly  produces  the  200 
kW(e)  that  is  desired.  The  thickness  of  the  electrical  power  plot  is  not  due  to  erratic 
behavior  of  the  electrical  circuit,  but  due  to  the  discrete  jumps  in  the  plot  from  the 
differences  in  the  cycle  averages.  MATLAB  connects  all  of  the  data  points  when 
plotting,  so  the  transition  between  cycles  produces  vertical  plot  lines  (this  is  shown 
clearly  in  Figure  6-11,  which  has  a  smaller  time  interval).  These  vertical  lines  can  be 
eliminated  with  lower  tolerances  for  the  numerical  solver,  but  some  error  was  introduced 
as  a  sacrifice  for  simulation  speed.  This  error  was  usually  less  than  two  percent  but 
occasionally  was  slightly  higher  than  three  percent. 

After  the  model  was  validated  at  steady  state,  perturbation  scenarios  were  created  for 
simulation  tests.  Table  6-5  lists  the  desired  perturbation  simulations  from  easiest  to 
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Figure  6-13.  System  performance  at  steady  state.  The  power  levels  are  approximately:  reactor 
(618  kW(t)),  Stirlings  (195  kW(e)),  radiators  (408  kW(t)). 

hardest.  The  scenarios  were  then  completed  in  this  order,  with  the  exception  of  the  last 
few  which  were  not  completed  due  to  time  constraints.  The  simulations  required  large 
amounts  of  time  when  they  were  run  on  a  personal  computer  (22X  simulation  time),  so 
all  system  simulations  were  run  on  a  multi-node  LINUX  cluster  in  the  Naval  Academy’s 
Computer  Aided  Design/Interactive  Graphics  (CADIG)  laboratory  (3X  simulation  time). 
The  rest  of  this  section  describes  the  results  of  these  simulations.  Due  to  the  large 
amount  of  data  available  from  each  simulation,  the  results  of  the  simulations  will  only 
show  the  reactor  power,  Stirling  electrical  power,  and  the  radiator  heat  rejection  so  that 
the  overall  system  performance  can  be  evaluated. 


Simulation 

System  Perturbations 

Variable  Altered 

0 

Steady  State 

none 

1 

0.001  Reactivity  Insertion 

reactivity 

2 

-0.001  Reactivity  Insertion 

reactivity 

3 

Radiator  efficiency  decrease 

background  temperature 

4 

Stirling  Failure 

number  of  converters 

5 

0.5  Ohm  Decrease  in  the  Resistive  Load 

load  resistance 

6 

1  Ohm  increase  in  the  Resistive  Load 

load  resistance 

7 

Load  short  out  (run  with  errors) 

load  resistance 

8 

Change  from  10%  to  100%  power  (was  not  completed) 

reactor  power 

9 

Change  from  100%  to  10%  power  (was  not  completed) 

reactor  power 

10 

System  start  up  (was  not  completed) 

reactor  power 

11 

System  shut  down  (was  not  completed) 

reactor  power 

Table  6-5.  List  of  planned  system  perturbation  simulations. 


System  at  Steady  State 


Reactor  Thermal  Power 
Stirling  Electrical  Power 
Radiator  Heat  Rejection 
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Simulation  1  was  a  0.001  reactivity  insertion  at  600  s  with  the  system  at  equilibrium. 

This  test  is  intended  to  show  what  the  effects  of  partially  removing  a  control  rod  are  on 
the  entire  system.  The  system  power  values  from  the  simulation  are  shown  in  Figure  6- 
14.  At  the  time  of  the  reactivity  insertion,  the  reactor  power  rapidly  increased  by  about 
20%.  As  a  result,  an  increased  thermal  load  was  added  the  Stirling  converters.  The 
Stirling  converters  responded  by  increasing  their  operating  power  and  the  radiators 


System  Response  to  a  0.001  Reactivity  Insertion  at  600  s 


Figure  6-14.  System  response  to  a  0.001  reactivity  insertion  at  600s.  The  final  power  levels  are 
approximately:  reactor  (640  kW(t)),  Stirlings  (205  kW(e)),  radiators  (420  kW(t)). 

rejected  more  heat.  The  Stirling  electrical  power  also  seemed  to  stabilize  at  about  610 
seconds.  This  is  due  to  the  fact  the  simulation  step  size  decreased  the  most  during  this 
sharp  change  in  the  transient.  This  result  showed  that  the  model  was  numerically  stable 
under  transient  conditions  and  predicts  that  such  a  system  would  be  load  following. 

Following  the  reactivity  insertion  simulation,  a  similar  simulation  was  run  with  a 
negative  0.001  reactivity  insertion  at  600  seconds  with  the  system  at  equilibrium.  The 


negative  reactivity  insertion  is  used  to  evaluate  how  partially  inserting  a  control  rod  in  the 


reactor  might  affect  the  overall  system.  This  simulation  (Simulation  2)  is  displayed  in 
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Figure  6-15.  Similar  to  the  previous  simulation,  the  reactor  responded  immediately. 
However,  in  this  case  the  reactor  power  rapidly  decreased  by  20%  due  to  the  negative 
reactivity.  The  resultant  decrease  in  thermal  load  on  the  Stirling  converters  caused  them 
to  operate  at  a  lower  power  level  (~  5%  decrease)  and  the  radiators  rejected  less  heat  (~ 
2.5%  decrease).  Again,  the  model  was  stable  and  predicted  that  the  system  would  be  load 
following. 


System  Response  to  a  -0.001  Reactivity  Insertion  at  600  s 


Figure  6-15.  System  response  to  a  negative  0.001  reactivity  insertion  at  600  s.  The  final  power  levels 
are  approximately:  reactor  (595  kW(t)),  Stirlings  (186  kW(e)),  radiators  (395  kW(t)). 

Simulation  3  was  a  radiator  efficiency  induced  perturbation.  This  perturbation 
simulates  a  loss  in  the  radiators  ability  to  reject  heat,  which  could  be  caused  by  a  change 
in  the  solar  or  albedo  incident  energy  on  the  radiation  panels.  Albedo  energy  is  solar 
energy  that  is  reflected  from  a  planet.  This  perturbation  was  modeled  by  increasing  the 
ambient  space  temperature  from  4  to  304  K  (see  Equation  5-5)  for  half  of  the  radiator 
panels,  which  creates  the  effects  of  approximately  4  kW/m  of  incident  energy  on  the 
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radiation  panels.  This  energy  flux  is  similar  to  the  solar  flux  the  panels  could  receive  if 
the  system  was  on  a  spacecraft  between  Mercury  and  Venus.  The  results  of  the 
simulation  are  shown  in  Figure  6-16.  The  inability  of  the  radiators  to  reject  heat  is  shown 
by  the  sudden  20%  drop  in  radiator  heat  rejection.  The  Stirling  converters  and  the  reactor 
both  respond  with  decreases  in  power  (~  3%  and  2%  decrease,  respectively),  which 
brings  the  system  back  into  equilibrium. 


System  Response  to  Change  in  Background  Temperature  from  4  to  304  K  at  600s 


Time  (s) 

Figure  6-16.  System  response  to  a  4  to  304  K  increase  in  ambient  space  temperature  for  one  side  of 
the  radiator  panels.  The  final  power  levels  are  approximately:  reactor  (596  kW(t)),  Stirlings  (186 

kW(e)),  radiators  (396  kW(t)). 

The  next  simulation  run  (Simulation  4)  was  a  Stirling  converter  failure.  Figure  6-17 
illustrates  how  the  reactor,  electrical,  and  radiated  power  of  the  system  responded  to  this 
perturbation.  As  the  converter  failed,  the  electrical  power  immediately  dropped  and  the 
entire  thermal  load  was  placed  on  the  three  operating  converters  instead  of  the  normal 
four.  The  three  remaining  converters  almost  accommodate  for  the  lost  converter  by 


increasing  power  levels.  The  reactor  and  radiators  responded  to  the  drop  in  power  by 


decreasing  power  levels.  As  time  increased  to  700  seconds,  the  reactor  power  partially 


83 


recovered  in  order  to  provide  more  heat  to  the  converters  operating  at  an  increased  power 
level.  Similarly  to  the  previous  simulations,  the  model  was  stable  and  predicted  that  the 
system  was  load  following. 


System  Response  to  a  Stirling  Engine  Failure  at  600  s 


Time  (s) 

Figure  6-17.  System  response  to  a  Stirling  converter  failure  at  600  s.  The  final  power  levels  are 
approximately:  reactor  (565  kW(t)),  Stirlings  (188  kW(e)),  radiators  (363  kW(t)). 

Next,  a  load  alteration  simulation  was  run  which  involved  cutting  the  resistive  load 
approximately  in  half  by  removing  half  an  ohm  of  resistance  (Rioad  in  Equation  4-65  was 
decreased  by  half  an  ohm).  In  this  simulation,  the  perturbation  simulates  the  effect  of  a 
significant  decrease  in  the  electrical  resistance  (possibly  through  short  or  turning 
instruments  on).  The  result  of  this  simulation  on  the  reactor  power,  radiator  heat 
rejection,  and  Stirling  electrical  power  is  shown  in  Figure  6-18.  The  decreased  resistance 
causes  an  increase  in  current,  which  increases  the  damping  on  the  piston.  The  increased 
damping  on  the  piston  caused  a  sudden  drop  in  Stirling  electrical  power,  but  the  higher 
current  draw  allowed  the  Stirling  converters  to  become  steady  at  200  kW(e)  while  using 
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less  thermal  energy.  In  this  case,  both  the  reactor  thermal  power  and  the  heat  rejected 
from  the  radiator  decreased  (both  ~  2%).  While  this  result  seems  incorrect,  it  can  be 
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Figure  6-18.  System  response  to  halving  the  load  resistance  at  600  s.  The  final  power  levels  are 
approximately:  reactor  (605  kW(t)),  Stirlings  (200  kW(e)),  radiators  (395  kW(t)). 


explained  by  the  changes  in  the  Stirling  converters’  heater  and  cooler  temperatures.  This 
caused  the  converters  to  slightly  increase  power  and  have  a  higher  efficiency  (~  34.7%). 
This  phenomenon  will  be  explained  further  in  the  explanation  of  the  next  simulation. 

The  last  simulation  completed  (Simulation  6)  involved  increasing  the  resistive  load  by 
adding  one  ohm  of  resistance  (Rioad  in  Equation  4-65  was  increased  by  one  ohm).  This 
simulation  was  meant  to  replicate  a  scenario  where  the  electrical  resistance  significantly 
increased  (due  to  bad  connections  or  suddenly  turning  off  components).  The  increased 
resistance  dramatically  lowered  the  current  in  the  electrical  circuit  and  consequently 
reduces  the  electrical  power  produced  by  the  system  (Figure  6-19).  The  sudden  loss  in 
current  led  to  a  sudden  loss  in  damping  which  allowed  the  piston  to  operate  at  a  higher 
speed.  However,  the  sudden  increase  (~  45%  of  original  value)  in  electrical  power  from 
this  perturbation  is  shortly  followed  by  a  decrease  (~  25%  of  original  value)  as  the 
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increase  in  heat  addition  from  the  Stirling  heater  to  the  working  gas  (Figure  6-20)  causes 
the  temperature  of  the  heater  to  drop  (Figure  6-21).  This  drop  in  heater  temperature 
eventually  causes  the  reactor  to  produce  more  power.  However,  as  shown  by  the  increase 
in  temperature  of  the  Stirling  cooler  (Figure  6-21),  the  Stirling  converters  are  now 
operating  less  efficiently  (~  21.7%).  In  others  words,  the  converters  are  requiring  more 
heat  input  to  produce  less  electrical  power,  meaning  the  system  is  operating  far  from  its 
optimum  efficiency  condition.  Simulations  5  and  6  demonstrate  that  the  electrical 
resistance  of  the  alternator  circuit  affects  the  efficiency  of  the  Stirling  converters.  Also, 
the  current  value  in  the  alternator  for  Simulation  6  experienced  a  shift  at  the  moment  of 
the  resistance  increase.  While  all  the  electrical  relationships  held  true  during  this  shift 
(including  electrical  power),  this  phenomenon  should  be  analyzed  and  better  understood. 


System  Response  to  1  Ohm  Resistance  Increase  at  600  s 


Time  (s) 

Figure  6-19.  System  response  to  doubling  the  load  resistance  at  600  s.  The  final  power  levels  are 
approximately:  reactor  (690  kW(t)),  Stirlings  (150  kW(e)),  radiators  (515  kW(t)). 
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Figure  6-20.  Stirling  heater  and  cooler  heat  transfer  values  during  a  doubling  of  the  electrical  load 
resistance.  The  heat  addition  of  the  heater  doubles  but  settles  at  a  value  13%  higher  than  the  initial 
value.  The  heat  rejection  of  the  cooler  increases  by  250%,  but  settles  at  a  value  20%  larger. 


Time  (s) 

Figure  6-21.  Stirling  heater  and  cooler  temperatures  during  a  doubling  of  the  electrical  load 
resistance.  The  heater  experiences  a  sudden  drop  in  temperature  while  the  cooler  experiences  a 

sudden  rise  in  temperature  at  600  s. 
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7.  Project  Conclusions 

At  the  conclusion  of  this  project,  a  model  of  a  Stirling  space  nuclear  power  system 
consisting  of  detailed  reactor,  Stirling,  and  radiator  submodels  was  developed. 

Simulations  of  the  individual  components  of  the  system  verified  that  the  assumptions 
used  were  reasonable  and  that  each  component  model  was  stable.  According  to  the 
simulations  presented  in  Section  6,  it  is  evident  that  the  system  modeled  tends  to  be  a 
stable,  load-following  system,  where  load  is  based  on  the  heat  transferred  and  rejected  by 
the  Stirling  heater  and  cooler  to  the  working  gas.  A  change  in  one  component  resulted  in 
changes  in  the  other  components  that  lead  the  system  to  reach  a  steady  state  condition.  It 
was  also  evident  that  each  part  of  the  system  had  a  different  time  lag  when  responding  to 
changes  in  the  system.  This  time  lag  is  inversely  proportional  to  the  components  mass 
and  size  (thermal  time  constant),  where  the  Stirling  converters  (smallest)  respond  the 
fastest  and  the  radiators  (largest)  respond  the  slowest.  For  Simulations  1-4,  the  power 
level  of  each  major  component  increased  or  decreased  in  the  same  direction  as  the  other 
components.  In  Simulations  5  and  6,  this  was  not  the  case.  For  the  unregulated  electrical 
system  modeled,  it  was  apparent  that  decreases  in  load  increased  thermal  efficiency  of  the 
Stirling  converters  while  increased  load  decreased  thermal  efficiency  (Figures  6-18  and 
6-19).  However,  in  all  cases,  the  reactor  power  matched  the  heat  drawn  by  the  Stirling 
converters  and  the  radiators  matched  the  heat  rejected  by  the  Stirling  converters. 

Conclusions 

•  Completed  model  was  stable 

•  Assumptions  used  were  reasonable 

•  Model  suggests  a  stable,  load  following  system 

•  Each  major  components  thermal  capacity  had  a  large  effect  on  the  systems 
dynamics 

•  Changing  electrical  load  affects  Stirling  converter  efficiency 
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This  project  has  much  potential  for  future  work.  First,  the  systems  individual 
elements  should  be  optimized  and  system  controls  should  be  examined  to  provide  a 
means  for  increasing  the  system’s  performance.  This  examination  should  include  the 
creation  of  a  detailed  performance  map  of  the  Stirling  converter  as  a  function  of  electrical 
load.  Second,  the  rest  of  the  simulations  listed  in  Table  6-5  should  be  run  and  further 
work  examining  the  effects  of  the  electrical  circuit  on  the  system  is  needed.  This  is 
especially  true  with  regard  to  the  current  anomaly  of  Simulation  6  in  Section  6.4.  Future 
work  should  also  study  the  effect  of  infinite  impedance  (opening  the  circuit).  A  test  with 
infinite  impedance  would  show  the  worst  case  scenario  for  the  circuit  resistance  since  it 
would  cause  the  largest  possible  increase  in  reactor  power.  Third,  a  more  detailed  and 
thorough  design  of  the  whole  system  would  provide  a  more  realistic  examination  of  the 
connections  between  the  three  main  components  and  the  performance/monetary  merits  of 
each  of  the  alternatives.  Fourth,  the  numerical  solvers  in  Simulink  should  be  examined 
further  so  that  a  thorough  analysis  of  accuracy  and  simulation  can  be  made.  Fifth,  the 
system  should  undergo  a  material  analysis  to  make  sure  that  all  of  the  components  can 
safely  withstand  operation  or  be  improved  with  the  selection  of  another  material.  Finally, 
it  is  recommended  that  the  model  be  compared  using  other  modeling  software  (such  as 
Simplorer®).  This  would  allow  for  comparisons  between  the  performance  of  each 
software  tool  and  it  could  serve  to  provide  an  assessment  as  to  the  adequacy  of  each 
model. 


Future  Work  Suggestions 

•  Further,  in-depth  sub-system  study  (including  a  Stirling  converter 
performance  map) 

•  Complete  the  planned  simulations  and  conduct  further  simulations  that 
examine  changes  in  the  alternator  electrical  circuit 

•  Develop  the  system  design  more  thoroughly,  especially  where  the  main 
components  connect 

•  Examine  Simulink’s  numerical  solvers  with  respect  to  accuracy  and  speed 

•  Conduct  a  material  evaluation  of  the  entire  system 

•  Compare  the  model  with  different  simulation  software 
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Appendix  A:  Reactivity  Insertion  of  0.001  at  the  5  Second  Mark.  Constant  600  kWt 

Heat  Flow  Out  of  the  Stirling  Hot  End. 
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of  the  reactor,  heat  transfer  loops,  and  heat  exchanger. 
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Appendix  B:  Reactivity  Insertion  of  -0.001  at  the  5  Second  mark.  Constant  600  kWt 

Heat  Flow  Out  of  the  Stirling  Hot  End. 

This  appendix  shows  all  the  output  values  for  the  title  simulation  from  the  partial  model 
of  the  reactor,  heat  transfer  loops,  and  heat  exchanger. 
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Appendix  C:  Ten  Degree  Temperature  Drop  for  Stirling  Hot  End  Heat  Exchanger. 

This  appendix  shows  all  the  output  values  for  the  title  simulation  from  the  partial  model 
of  the  reactor,  heat  transfer  loops,  and  heat  exchanger. 
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Appendix  D:  Ten  Degree  Temperature  Rise  for  Stirling  Hot  End  Heat  Exchanger. 

This  appendix  shows  all  the  output  values  for  the  title  simulation  from  the  partial  model 
of  the  reactor,  heat  transfer  loops,  and  heat  exchanger. 
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Appendix  E:  Stirling  Model  Data  Provided  by  NASA 


TO:  Justin  Langlois 

FROM:  Barry  Penswick 

SUBJECT:  50  KW  Stirling  convertor  model 

DATE:  6  Dec  05 

The  following  material  describes  the  parameters  of  a  hypothetical  50  KW  Stirling  cycle 
convertor  that  hopefully  will  meet  your  system  modeling  requirements.  Some  of  the 
generic  convertor  characteristics  can  be  related  back  to  a  relatively  recent  review  carried 
out  at  NASA  /  GRC  concerning  high  power  space  systems.  As  such  the  operating  point 
employed  in  the  model  may  not  match  potential  operating  conditions  that  may  evolve 
from  your  ongoing  work,  however  modifying  the  model  to  meet  your  specific 
requirements  is  a  straightforward  process. 

Configuration  Issues 

The  reference  configuration  of  the  50  KW  system  is  a  conventional  piston  /  displacer 
arrangement  with  a  linear  alternator  attached  to  the  power  piston.  The  basic  external 
configuration  of  the  convertor  is  shown  in  the  following  sketch.  It  is  important  to  note  that 
as  configured  the  50  KW  convertor  is  NOT  balanced.  I’ve  assumed  that  you  will  have  to 
address  this  issue  and  would  recommend  that  you  simply  assume  that  the  two 
convertors  are  operated  in  a  “hot  end  to  hot  end”  opposed  configuration  and  are 
electrically  coupled  so  as  to  provide  the  proper  synchronization  of  the  pistons  and  in  turn 
effective  vibration  balancing  the  pair  of  50  KW  convertors. 


Rejector  HX  Assembly 


Heater  Assembly 


Regenerator 


Alternator  Housing  &  PV 


Figure  1.  Convertor  configuration 
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From  the  viewpoint  of  the  internal  component  configuration  it  is  assumed  that  the 
conventional  common  cylinder  diameter  with  the  displacer  rod  through  the  power  piston 
is  employed.  The  regenerator  is  of  an  annular  configuration  surrounding  the  displacer. 
Both  the  accepter  and  rejecter  heat  exchangers  are  of  the  tubular  type.  This 
fundamental  configuration  is  currently  employed  by  SunPower  Inc.  and  other  groups  and 
represents  a  logical  configuration  for  this  sized  convertor.  Terminology  employed  in  the 
model  is  noted  on  the  following  sketch  where  applicable  and  values  provided  in  the  next 
table. 


Figure  2.  Convertor  internal  component  configuration  and  terminology 
Table  1  -  Convertor  parameters 


Parameter 

Nominal  Value  /  Comment 

General 

Operating  frequency 

90  Hz 

Charge  pressure 

1 36  Bar 

Effective  hot  end  temperature 

925  K  (652  C) 

Effective  rejecter  temperature 

463  K  (190  C) 

Total  convertor  mass 

270  Kg 

Net  electrical  out  put  power 

50  KW 

Convertor  cvcle  efficiency 

.316 

Convertor  over  all  efficiency 

.294 

Power  piston  related 

Piston  frontal  area 

.0381  mA2 

Piston  amplitude 

13  mm 

Power  piston  /  alternator  moving  mass 

25  Kg 
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Compression  space  volume 

6.0E-4  mA3 

Bounce  space  volume 

14.5E-3  mA3 

Piston  /  displacer  phase  angle 

59.3  degrees 

Piston  PV  power 

53.6  KW 

Displacer  related 

Displacer  hot  end  frontal  area 

3.892E-2 

Displacer  cold  end  frontal  area 

3.81  E-2 

Displacer  rod  diameter 

32  mm 

Displacer  amplitude 

8.4  mm 

Displacer  effective  moving  mass 

8.5  Kg 

Expansion  space  volume 

4.34E-4  mA3 

Rejecter  HX  related 

Length 

75  mm 

Tube  ID 

1 .4  mm 

Tube  number 

3100 

Regenerator  related 

Frontal  area 

5.01 7E-2 

Length 

55  mm 

Porosity 

88% 

Wire  diameter 

25  micron 

Accepter  HX  related 

Length 

75  mm 

Tube  ID 

1 .3  mm 

Tube  number 

2750 

Linear  Alternator 


A  simplified  alternator  modeling  technique  was  employed  to  develop  a  reference 
alternator.  Since  there  was  no  definition  of  the  alternator  requirements  the  values  noted 
in  the  following  table  represent  the  parameters  of  the  reference  alternator.  As  noted 
previously  it  is  a  relatively  straightforward  process  to  adapt  the  model  alternator  to  your 
specific  operating  requirements. 

Table  2  -  Linear  alternator  parameters 


Alternator  Characteristics 

Nominal  Value  /  Comments 

Design  output  power 

50  KW 

Terminal  voltage 

200  Vrms 

Design  efficiency 

93% 

Operating  temperature  (apx.) 

490  K  (215  C) 

Piston  amplitude 

13  mm 

Alternator  constant 

38.5  V  sec/m 

Coil  resistance  (est.) 

.028  ohm 

Coil  inductance  (est.) 

.2  mH 

Appendix  F :  MATLAB®  Code  for  System  Constants  and  Initial  Conditions 


%  Fuel  Properties  (UN) 

beta=0.0065;  %delayed  neutron  fraction  (dl) 
lamda=0.077;  %one  group  precursor  decay  constant  (1/s) 
ng=10A(-6);  %neutron  generation  time  (s) 
rowo=0;  %initial  reactivity  (dl) 

Lf=0.25;  %fuel  pin  length  (m) 
kf=14;  %fuel  conductivity  UN  (W/ (m*K)  ) 
mU235=180;  %U-235  mass  (kg) 
mUN=2*mU235* ( 235+15 ) /235 ;  %UN  mass  (kg) 
pf=14320;  %fuel  density  UN  (kg/mA3) 

N=108;  %number  of  fuel  pins  (dl) 
vf=mUN/pf;  %fuel  volume  (mA3 ) 

rf =round ( 10000*sqrt ( vf / (pi*Lf *N) ) ) /10000 ;  %fuel  pin  radius  (m) 
cpf=219.4;  %fuel  specific  heat  UN  (W*s/ (kg*K) ) 
alphaf=-3*10A (-5) ;  %fuel  temperature  coefficient  U02  (1/K) 

%  Power  level 

Po=6.9*10A5;  %initial  power  level  (Wt) 

PEo=Po*beta/ ( lamda*ng) ;  %initial  power  equivalence  (Wt) 

%  Moderator  Properties  (NaK) 

pm=725;  %moderator  density  (kg/mA3) 

cpml=875;  %mod.  spec,  heat  at  reactor  inlet  (W*s/ (kg*K) ) 

cpm2=875;  %mod.  spec,  heat  at  reactor  outlet  (W*s/ (kg*K) ) 

cpm= ( cpml+cpm2 ) /2 ;  %avg.  mod.  spec,  heat  (W*s/ (kg*K) ) 
mdotm=Po/ ( 2*cpm*25 ) ;  %mod.  mass  flow  (kg/s) 
pc=sqrt (2*pi) *rf ;  %pitch  (m) 

Acell=pcA2-pi*rf A2 ;  %cell  mod.  area  (mA2) 
pw=pi*2*rf;  %wetted  perimeter  (m) 
vm=Acell*Lf *N;  %moderator  volume  (mA3) 
de=4*Acell/pw;  %equivalent  diameter  (m) 
alpham=0;  %mod.  temp,  coeff.  (1/K) 
velm=mdotm/ (pm*Acell*N) ;  %mod.  velocity  (m/s) 
um=l . 59*10A (-4) ;  %mod.  viscosity  (kg/(m*s)) 
km=25.6;  %mod.  thermal  conductivity  (W/ (m*K) ) 

Re=de*velm*pm/um;  %mod.  Reynolds  number  (dl) 

Pr=cpm*um/km;  %mod.  Prandtl  number  (dl) 

Pe=Re*Pr;  %Peclet  Number  (dl) 

ht= (km/de) * ( 7+4 . 24* (pc/ (2*rf ) ) A1 . 52) ;  %mod.  heat  transfer  coeff.  (W/ (mA2*K) ) 

%  Primary  Loop  Outside  the  Reactor 
dp=0.04;  %pipe  diameter  (m) 

Ap=pi*dpA2/4;  %pipe  cross-section  area  (mA2) 

tau3=2.5;  %time  delay  from  reactor  outlet  to  heat  exchanger  inlet  (s) 
tau4=2.5;  %time  delay  from  heat  exchanger  outlet  to  reactor  inlet  (s) 
lg3=mdotm*tau3/ (pm*Ap) ;  %pipe  length  from  reactor  outlet  to  hxl  (m) 
lg4=mdotm*tau4/ (pm*Ap) ;  %pipe  length  from  hxl  to  reactor  inlet  (m) 

%  Primary  to  Secondary  Heat  Exchanger  Properties  (steel) 

dt=0.0095;  %diameter  of  hx  tubes 

tt=0.003;  %hx  tube  thickness 

odt=dt+2*tt;  %outer  diameter  tube 

At=pi*dt A2/4;  %area  cs  of  hx  tubes 

n=18;  %hx  size 

Nt=nA2+4*n;  %number  of  hx  tubes 
Lt=0.5;  %tube  length 

vhxl=At*Nt*Lt ;  %volume  of  coolant  in  heat  exchanger  (mA3) 
ahxl=pi*dt*Lt*Nt ;  %primary  side  heat  exchanger  wall  area  (mA2) 

Rem=mdotm*dt/ (At*Nt*um) ;  %hxl  Reynolds  number 
Prm=cpm*um/km;  %hxl  Prandtl  number 
Pem=Rem*Prm;  %hxl  Peclet  number 

hxl= (km/dt) * ( 4 . 82+0 . 697*odt/dt ) ;  %primary  side  heat  transfer  coeff.  (W/ (mA2*K) ) 

ahx2=Lt*pi*odt*Nt ;  %wall  area  loop  2 

phx2=l . 25*odt ;  %hx  pitch 

pw=8000;  %wall  density  (kg/mA3) 

vw=pi/4* (odtA2-dtA2 ) *Lt*Nt ;  %wall  volume  (mA3) 

mw=vw*pw;  %wall  mass  (kg) 

cpw=500;  %wall  spec,  heat  steel  (W*s/ (kg*K) ) 
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kw=21.5;  %wall  conductivity 
Bixl=hxl*Lt/kw;  %tube  Biot  number 

mdots=mdotm;  %mass  flow  rate  of  coolant  in  secondary  loop  (kg/s) 

ps=pm;  %density  of  sec.  coolant  (kg/mA3) 

cps=cpm;  %avg.  spec,  heat  of  sec.  coolant  (W*s/ (kg*K) ) 

cpsl=cpm;  %avg.  spec,  heat  of  sec.  coolant  at  exch.  inlet  (W*s/(kg*K)) 
cps2=cpm;  %avg.  spec,  heat  of  sec.  coolant  at  exch.  outlet  (W*s/(kg*K)) 
ks=km;  %sec.  coolant  conductivity 
us=um;  %sec.  coolant  viscosity 
ds=2*n*phx2;  %  shell  diameter  (m) 

As=pi/4*dsA2-Nt*pi/4*odt A2 ;  %shell  area  (mA2) 
dsh=4*As/ (Nt*pi*odt ) ;  %hydraulic  diameter  of  shell 
Res=mdots*ds/ (ahx2*us) ;  %hx2  Reynolds  number 
Prs=cps*us/ks ;  %hx2  Prandtl  number 
Pes=Res*Prs;  %hx2  Peclet  number 

hx2= (ks/dsh) * ( 7+4 . 24* (phx2/odt) A1 . 52) ;  %secondary  side  heat  tranfer  coeff.  (W/ (mA2*K) ) 
Bix2=hx2*Lt/kw;  %shell  Biot  number 

vhx2=Lt*As;  %volume  of  sec.  coolant  in  heat  exch.  (mA3) 

UAhx=l/ (1/ (hxl*ahxl) +log(odt/dt) / (2*pi*kw*Lt*Nt) +1/ (hx2*ahx2) ) ;  %hx  total  resistance 
dThx=Po/UAhx;  %temperature  drop  across  heat  exch.  (K) 

%  Secondary  Loop 

tau6=2.5;  %time  delay  from  heat  exchanger  outlet  to  Stirling  engine  heat  source  (s) 
tau5=2.5;  %time  delay  from  Stirling  heat  sink  to  heat  exchanger  inlet  (s) 

%  Stirling  Heat  Source 

%Heat  Pipe  -  CP  Grade  2  Titanium 

HHxLg=0.075;  %Length  of  the  Stirling  Hot  Heat  Exchanger  (m) 

HHxA=5 . 017*10A (-2 ) ;  %Hot  Heat  Exchanger  Cross  Section  Area  (mA2) 

nht=2750*4;  %number  of  hot  heat  exchanger  gas  tubes 

htd=0 . 0013/2 ;  %Hot  Heat  Exchanger  gas  tube  diameter  (m) 

htA=nht*pi*htdA2/4 ;  %Hot  Heat  Exchanger  gas  tubes  area  (mA2) 

sf=0.5;  %heat  pipe  area  safety  factor 

hpid=sf *0 . 05 ;  %heat  pipe  inner  diameter  (m) 

hpt=sf *7 . 11*10A (-4) ;  %heat  pipe  thickness  (m) 

hpod=hpid+2*hpt ;  %heat  pipe  outer  diameter  (m) 

hpwt=sf *0 . 0015 ;  %heat  pipe  wick  thickness  (m) 

hpwid=hpid-2*hpwt ;  %heat  pipe  wick  inner  diameter  (m) 

hpposs= (HHxA-htA) / (pi*hpodA2/4 ) ;  %number  of  heat  pipes  possible 

Le=l;  %heat  pipe  evaporator  length  (m) 

Lc=HHxLg;  %heat  pipe  condenser  lenght  (m) 

kp=22;  %heat  pipe  conductivity  -  titanium  (W/ (mA2*K) ) 

Cp=540;  %heat  pipe  specific  heat  -  titanium  ( J/ (kg*K) ) 
rowp=4510;  %heat  pipe  density  -  titanium  (kg/mA2) 
porp=0.65;  %heat  pipe  porosity 

kl=32.3;  %heat  pipe  fluid  conductivity  -  potassium  (W/ (m*K)  ) 

kef  f =kl*  (  (kl+kp)  -  ( 1-porp)  *  (kl-kp)  )  /  (  (kl+kp)  +  ( 1-porp)  *  (kl-kp)  )  ;  %heat  pipe  wick  effective 
thermal  conductivity  (W/ (m*K) ) 

Rpe=log (hpod/hpid) / (2*pi*Le*kp) ;  %heat  pipe  evaporator  thermal  resistance  (K/W) 

Rwe=log (hpid/hpwid) / ( 2*pi*Le*kef f ) ;  %heat  pipe  evaporator  wick  thermal  resistanc  (K/W) 
Rpc=log (hpod/hpid) / ( 2*pi*Lc*kp) ;  %heat  pipe  condenser  thermal  resistance  (K/W) 

Rwc=log (hpid/hpwid) /( 2*pi*Lc*kef f ) ;  %heat  pipe  condenser  thermal  resistance  (K/W) 
Rtot=Rpe+Rwe+Rpc+Rwc;  %heat  pipe  total  thermal  resistance  (K/W) 

Nhp=f loor (hpposs ) ;  %number  of  heat  pipes 
NSE=4;  %number  of  Stirling  engines 

Apxs=pi* (hpodA2-hpidA2 ) /4 ;  %heat  pipe  cross  section  area  (mA2) 

Vpe=Apxs*Le*Nhp*NSE;  %heat  pipe  condenser  volume  (mA3) 

Aevap=pi*hpod*Le*Nhp*NSE;  %heat  pipe  evaporator  transfer  area  (mA2) 
hppitch=2*hpod;  %heat  pipe  matrix  pitch  (m) 

hpf lowarea= (hppitchA2-pi*hpodA2/4 ) ;  %heat  pipe  unit  cell  flow  area  (mA2) 
hpwetper=pi*hpod;  %heat  pipe  unit  cell  wetted  perimeter  (m) 
ahx3=hpf lowarea*Nhp;  %total  flow  area  (mA2) 
vhx3=ahx3*Le;  %flow  volume  (mA3) 

hphyddiam=4*hpf lowarea/hpwetper ;  %heat  pipe  unit  cell  hydraulic  diameter  (m) 
uf low=mdots/ (hpf lowarea*ps*4 ) ;  %flow  velocity  (m/s) 

Red=ps*uf low*hphyddiam/um;  %flow  Reynolds  number 
Nux3=7+4 . 24* ( 2 ) A1 . 52 ;  %flow  Nusselt  number 

hx3=Nux3*ks/hphyddiam;  %heat  xfer  coef.  between  flow  and  heat  pipes  (W/ (mA2*K) ) 
deltaT=Po*Rtot/ (Nhp*NSE) ;  %temperature  difference  between  flow  and  Stirling  hot  heat 
exchanger  (K) 

mh=20;  %heater  mass  -  TeCu  (kg) 
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cph=393.5;  %heater  specific  heat  -  TeCu  ( J/ (kg*K) ) 

%Temperatures  (K) 

Twh=925;  %Stirling  heater 
Tpe=Twh+deltaT;  %heat  pipe  evaporator 

Td2o=Tpe+ (Po/ (mdots*cps) ) / (exp (hx3*Aevap/ (mdots*cps) ) -1 ) ;  %init.  coolant  temp,  at  heat 
source  outlet  (K) 

Td2=Tpe+ (Po/ (mdots*cps) ) / (exp (hx3*Aevap/ (mdots*cps) ) -1) ;  %coolant  temp,  at  heat  source 
outlet  (K) 

Tdlo=Td2o+Po/ (mdots*cps) ;  %init.  avg.  coolant  temp,  at  heat  source  (K) 

Tdl=Td2+Po/ (mdots*cps) ;  %avg.  coolant  temp,  at  heat  source  (K) 

Tclo=Td2;  %init .  temp,  of  sec.  fluid  entering  exch.  (K) 

Tcl=Td2;  %temp.  of  sec.  fluid  entering  exch.  (K) 

Tco= (Tdlo+Tclo) /2 ;  %initial  avg.  temp,  of  secondary  fluid  at  exch. (K) 

Tc= (Tdl+Tcl ) /2 ;  %avg.  temp,  of  secondary  fluid  at  exch. (K) 

Tc2=Tc+Po/ ( 2*mdots*cps ) ;  %temp.  of  sec.  fluid  leaving  exch.  (K) 

Tbo=Tco+Po/UAhx;  %init .  avg.  temp,  of  heat  exchanger  primary  coolant  (K) 

Tb=Tc+Po/UAhx;  %avg.  temp,  of  heat  exchanger  primary  coolant  (K) 

Tb2o=Tbo-Po/ ( 2*mdotm*cpm) ;  %init.  heat  exchanger  primary  coolant  outlet  temp.  (K) 
Tb2=Tb-Po/ ( 2*mdotm*cpm) ;  %heat  exchanger  primary  coolant  outlet  temp.  (K) 

Tsl=Tbo- ( 1/ (hxl*ahxl ) ) *UAhx*dThx;  %initial  wall  temp.  (K) 

Ts2=Tco+ ( 1/ (hx2*ahx2 ) ) *UAhx*dThx;  %wall  temp.  (K) 

Tblo=2*Tbo-Tb2 ; %init .  coolant  temp,  at  heat  exchanger  inlet  (K) 

Tbl=2*Tb-Tb2 ; %coolant  temp,  at  heat  exchanger  inlet  (K) 

Tmlo=2*Tbo-Tblo; %init .  mod.  temp,  at  reactor  outlet  (K) 

Tml=2*Tb-Tbl ;  %mod.  temp,  at  reactor  outlet  (K) 

Tmo=Tmlo+Po/ ( 2*mdotm*cpm) ; %initial  avg.  mod.  temp.(K) 

Tm=Tml+Po/ ( 2*mdotm*cpm) ;  %mod.  avg.  temp.  (K) 

Tm2=2*Tm-Tml ;  %reactor  outlet  coolant  temp.  (K) 

Tf o=Tmo+ ( l+2*kf / (rf *ht ) ) *Po/ ( 4*pi*Lf *kf *N) ; %initial  fuel  temperature  (K) 

Tf=Tm+ ( l+2*kf / (rf *ht ) ) *Po/ ( 4*pi*Lf *kf *N) ; %fuel  temperature  (K) 

%Stirling  Parameters  and  IC's 

vc=1.3;  %compression  space  volume  factor  -  multiple  of  volume  given  by  NASA 
ve=l;  %expansion  space  volume  factor  -  multiple  of  volume  given  by  NASA 
Xp=0.013;  %piston  displacement  amplitude  (m) 

Xd=0.0084;  %displacer  displacement  amplitude  (m) 
mp=20.3;  %piston  mass  (kg) 
md=8.5;  %displacer  mass  (kg) 

Dd=0;  %displacer  damping  (kg/s) 

kdisp=3 . 1*10 A6 ;  %displacer  spring  constant  (N/m) 
kd3=0*10A10;  %third  order  displacer  spring  constant  (N/mA3) 
kd=kdisp-kd3*XdA2 ;  %ad justed  first  order  spring  constant  (N/m) 
kalt=38.5;  %alternator  constant  (V*s/m) 

f req=sqrt (kd/md) / ( 2*pi ) ;  %expected  operating  frequency  (Hz) 

Vmax=kalt*Xp* ( 2*pi*f req) ;  %max  voltage  expected  (V) 
kWe=50000;  %design  power  output  per  Stirling  (kWe) 

Zalt=VmaxA2/ (2*kWe) ;  %alternator  impedance  (ohms) 

Apist=0 . 0381 ;  %piston  area  (mA2) 
xpo=0;  %initial  piston  position  (m) 
vpo=Xp*2*pi*f req;  %initial  piston  velocity  (m/s) 
kalt=38.5;  %alternator  constant  (Vs/m) 

Ralt=0.028;  %alternator  resistance  (ohm) 

Halt=0.0002;  %alternator  inductance  (H) 

Calt=0;  %circuit  capacitance  (F)  -  power  corrected  value  =  1/ ( ( 2*pi*f req) A2*Halt ) ; 
Rload=sqrt (ZaltA2- (2*pi*freq*Halt) A2) -Ralt;  %electric  load  resistance  (ohms) 
Ialt=kalt*vpo/Zalt ;  %max  expected  current  (A) 

Dp=kalt A2/Zalt ;  %piston  damping  (kg/s) 
drd=0.032;  %displacer  rod  diameter  (m) 

Arod=pi*drdA2/4 ;  %displacer  rod  area  (mA2) 

Ade=0. 03892;  %displacer  area-exp.  side  (mA2) 

Adc=Apist;  %displacer  area-comp,  side  (mA2) 

dpist=sqrt ( 4*Ade/pi) ;  %piston  diameter  (m) 

phase=60;  %expected  piston-displacer  phase  angle  (deg) 

xdo=Xd*sin (pi*phase/180 ) ;  %displacer  initial  position  (m) 

vdo=Xd*2*pi*freq*cos (pi*phase/180 ) ;  %displacer  initial  velocity  (m/s) 

pmean=13 . 6*10A6 ;  %piston  mean  pressure  (Pa) 

Twk=463;  %cooler  temperature  (K) 

Twh=925;  %heater  temperature  (K) 

Twr= (Twh+Twk) /2 ;  %regenerator  temp.  (K) 

Vcso=vc*6*10 A (-4 ) ;  %compression  space  volume  (mA3) 
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Nkchan=3100*4 ;  %number  cold  exchanger  channels  (dl) 

CHxLg=0.075;  %cold  exch.  length  (m) 
dtk=0 . 0014/2 ;  %cooler  channel  diameter  (m) 

CHxpw= (pi*dtk) *Nkchan;  %cooler  wetted  perimeter  (m) 

Awk=CHxpw*CHxLg;  %cooler  wall  area  (mA2) 

ACHx=Nkchan*pi*dtkA2/4 ;  %cooler  cross-sectional  area  (mA2) 

CHxDh=4*ACHx/CHxpw;  %cooler  hydraulic  diameter  (m) 

Vk=ACHx*CHxLg;  %cooler  volume  (mA3) 
rLg=0.055;  %regen.  1  length  (m) 
por=0.88;  %regen.  porosity  (dl) 

Acr=5 . 017*10A (-2 ) ;  %regen.  area  (mA2) 

Ar=Acr*por;  %regen.  flow  area  (mA2) 

dr=2/pi* (pi* (Apist+Acr ) ) A ( 1/2 ) ;  %outer  regen.  diameter  (m) 
dir=sqrt ( 4*Ade/pi ) ;  %inner  regen.  diameter  (m) 

Awr=pi* (dr+dir ) *rLg;  %regen.  1  wall  area  (mA2) 
rwd=25*10 A (-6 ) ;  %regen.  wire  diameter  (m) 
rpw=4* ( 1-por ) *Acr/rwd;  %reg.  wetted  perimeter  (m) 
rDh=4*Ar /rpw;  %regen.  hyd.  diameter  (m) 

Amr=rpw*rLg;  %regen.  1  matrix  area  (mA2) 

Vr=rLg*Ar;  %regen.  1  volume  (mA3) 

HHxLg=0.075;  %heater  length  (m) 
htd=. 0013/2;  %heater  tube  diameter  (m) 

Nht=2750*4;  %number  heater  tubes  (dl) 

AHHx=Nht*pi/4*htdA2 ;  %heater  cross-sectional  area  (mA2) 

Vh=AHHx*HHxLg;  %heater  volume  (mA3) 

HHxpw=Nht*pi*htd;  %heater  wetted  perimeter  (m) 

HHxDh=4*AHHx/HHxpw;  %heater  diameter  hyd.  (m) 

Awh=HHxpw*HHxLg;  %heater  wall  area  (mA2) 

Veo=ve*4 . 34*10 A (-4 ) ;  %expansion  space  volume  (mA3) 

Vbso=14 . 5*10 A (-3 ) ;  %bounce  space  volume  (mA3) 

Vt=Vcso+Veo+Vh+Vr+Vk;  %total  volume  (mA3) 
gammaHe=l . 66 7;  %helium  ratio  of  specific  heats  (dl) 

Rgas=2077;  %helium  gas  constant  ( J/ (kg*K) ) 
kHe=0.14;  %helium  conductivity  ( J/ (kg*K) ) 
cp=5230;  %helium  specific  heat-const,  press.  ( J/ (kg*K) ) 
cv=cp-Rgas;  %helium  specific  heat-const,  vol .  ( J/ (kg*K) ) 

cs=500;  %steel  specific  heat  ( J/ (kg*K) ) 

CSutherland=l . 286e-6 ;  %Sutherland ' s  law  constant 

SSutherland=8 . 63 ;  %Sutherland ' s  law  constant 

Diss=0;  %initial  total  Stirling  gas  energy  dissipation  (W) 

Dissk=0;  %initial  Stirling  cooler  gas  energy  dissipation  (W) 

Dissr=0;  %initial  Stirling  regen.  gas  energy  dissipation  (W) 

Dissh=0;  %initial  Stirling  heater  gas  energy  dissipation  (W) 

pcso=l . 3243e+007;  %initial  compression  space  pressure  (Pa) 

peo=l . 3174e+007;  %initial  expansion  space  pressure  (Pa) 

pko=pmean;  %initial  cooler  pressure  (Pa) 

pro=pmean;  %initial  regen.  pressure  (Pa) 

pho=pmean;  %initial  heater  pressure  (Pa) 

WpV=53 . 6*10 A3 ;  %expected  Stirling  cycle  PV  power  (W) 
cycef f=0 . 316 ;  %expected  cycle  efficiency 

dQhss=WpV/cycef f ;  %expected  cycle  avg.  heat  addition  (W) 
dQkss=WpV-dQhss;  %expected  cycle  avg.  heat  rejection  (W) 

Nuh=3.66;  %initial  Stirling  heater  Nusselt  number 
Nuk=3.66;  %initial  Stirling  cooler  Nusselt  number 

hAk=Nuk*kHe*Awk/CHxDh;  %initial  Stirling  cooler  heat  transfer  coefficient  (W/ (mA2*K) ) 
hAh=Nuh*kHe*Awh/HHxDh;  %initial  Stirling  heater  heat  transfer  coefficient  (W/ (mA2*K) ) 
hAr=por A2 . 61*kHe*Awr/rDh;  %initial  Stirling  regen.  heat  transfer  coefficient  (W/ (mA2*K) ) 
Th=Twh-dQhss/hAh;  %initial  heater  space  temp.  (K) 

Tk=Twk-dQkss/hAk;  %initial  cooler  space  temp.  (K) 

Tcs=Tk;  %initial  compression  space  temperature  (K) 

Te=Th;  %initial  expansion  space  temperature  (K) 

Tr= (Th-Tk) /log (Th/Tk) ;  %initial  regen.  temp.  (K) 

Tmr=Tr;  %initial  regen.  matrix  temp.  (K) 

muk=CSutherland*TkAl . 5/ (Tk+SSutherland) ;  %initial  cooler  viscosity 
muh=CSutherland*ThAl . 5/ (Th+SSutherland) ;  %initial  heater  viscosity 
mur=CSutherland*Tr A1 . 5/ (Tr+SSutherland) ;  %initial  regenerator  viscosity 
meo=peo* (Veo-xdo*Ade) / (Rgas*Twh) ;  %initial  expansion  space  gas  mass  (kg) 
mho=pho*Vh/ (Rgas*Th) ;  %initial  heater  space  gas  mass  (kg) 

mro=pro*Vr/ (Rgas*Tr ) ;  %initial  regen.  space  gas  mass  (kg) 

mko=pko*Vk/ (Rgas*Tk) ;  %initial  cooler  space  gas  mass  (kg) 

mcso=pcso* (Vcso+ (xdo-xpo) *Apist ) / (Rgas*Twk) ;  %initial  compression  space  gas  mass  (kg) 
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mbo=pmean* (Vbso+xdo*Arod+xpo*Apist ) / (Rgas*Twk) ;  %initial  bounce  space  gas  mass  (kg) 

mwork= (meo+mho+mro+mko+mcso) ;  %total  gas  mass  of  the  working  spaces  (kg) 

dQh=hAh* (Twh-Th) ;  %initial  heater  heat  transfer  (W) 

dQk=hAk* (Twk-Tk) ;  %initial  cooler  heat  transfer  (W) 

dQr=hAr* (Twr-Tr ) ;  %initial  regen.  heat  transfer  (W) 

dQ=dQh+dQr+dQk;  %initial  total  heat  transfer  (W) 

%Radiator  Panels 

Chanht=0 . 0725;  %cooler  NaK  channel  height  (m) 

Chanwd=0.1;  %cooler  NaK  channel  width  (m) 

ChanLg=pi*dpist+2* (dr-dpist ) ;  %cooler  NaK  channel  length  (m) 

NChan=2;  %number  of  cooler  NaK  channels 
ptl=840;  %NaK  density  in  rejection  loop  (kg/mA3) 
cptl=915;  %NaK  specific  heat  in  rejection  loop  ( J/ (kg*K) ) 
ktl=25;  %NaK  conductivity  in  rejection  loop  (W/ (m*K)  ) 

AChan=Chanht*Chanwd;  %NaK  channel  cross  sectional  area  (mA2) 

VChan=NChan*AChan*ChanLg;  %NaK  channel  volume  (mA3) 

PChan=2*Chanht+2*Chanwd;  %NaK  channel  wetted  perimeter  (m) 

ChanDh=4*AChan/PChan;  %NaK  channel  hydraulic  diameter  (m) 

AwChan=NChan*PChan*ChanLg;  %NaK  channel  wetted  area  (mA2) 
mdotChan=mdotm/4 ;  %NaK  channel  mass  flow  (kg/s) 
mdott=mdotm;  %rejection  loop  mass  flow  (kg/s) 

Arlp=NChan*AChan*4 ;  %rejection  loop  pipe  cross  sectional  area  (mA2) 
drlp=sqrt ( 4* (Arlp) /pi) ;  %rejection  loop  pipe  diameter  (m) 
dradhp=0 . 0254 ;  %radiator  heat  pipe  diameter  (m) 

idradhp=dradhp-2*0 . 049*0 . 0254;  %radiator  heat  pipe  inner  diameter  (m) 
idhpw=idradhp-2*0 . 049*0 . 0254 ;  %radiator  heat  pipe  wick  inner  diameter  (m) 
utl=3 . 81*10 A (-4 ) ;  %NaK  viscosity  in  rejection  loop  (kg/(m*s)) 

PrChan=utl*cptl/ktl;  %NaK  channel  Prandtl  number 

velChan=mdotChan/ (ptl*AChan*NChan) ;  %NaK  channel  flow  velocity  (m/s) 
veltl=mdott/ (ptl*Arlp) ;  %rejection  loop  flow  velocity  (m/s) 

ReChan=ptl*velChan*ChanDh/utl;  %NaK  channel  Reynolds  number 
PeChan=PrChan*ReChan;  %NaK  channel  Peclet  number 

hChan=ktl/ChanDh* ( 9 . 49+0 . 0596*PeChanA0 . 86 ) ;  %NaK  channel  heat  transfer  coefficient 
(W/ (mA2*K) ) 

dQChan=4 . 06e+005/4 ;  %NaK  channel  heat  transfer  (W) 

dTChan=dQChan/ (hChan*AwChan) ;  %temperature  difference  between  cooler  and  NaK  channel  (K) 
dTax=dQChan/ (mdotChan*cptl ) ;  %temperature  drop  across  NaK  channel  (K) 
Rehp=ptl*veltl*dradhp/utl;  %flow  Reynolds  number  at  radiator  heat  pipes 
Pehp=PrChan*Rehp;  %flow  Peclet  number  at  radiator  heat  pipes 

hradhp= (ktl/dradhp) * (5+0 . 025*PehpA0 . 8) ;  %heat  transfer  coefficient  at  radiator  heat  pipes 
(W/ (mA2*K) ) 

Lerhp=drlp*2 ;  %radiator  heat  pipe  evaporator  length  (m) 
dwrhp=0 . 0027*0 . 0254;  %radiator  heat  pipe  wetted  diameter  (m) 
erhp=0.2738;  %radiator  heat  pipe  wick  porosity 
kwhp=0.678;  %radiator  heat  pipe  water  conductivity  (W/ (m*K) ) 

kef f rhp=kl* ( (kwhp+kp) - ( 1-erhp) * (kwhp-kp) ) / ( (kwhp+kp) + ( 1-erhp) * (kwhp-kp) ) ;  %radiator  heat 
pipe  effective  wick  conductivity  (W/ (m*K) ) 

Aerhp=Lerhp*dradhp*pi;  %radiator  heat  pipe  evaporator  area  (mA2) 
sigmasb=5 . 67*10 A (-8 ) ;  %Stephan-Boltzman  constant  (W/ (mA2*KA4) ) 
emiss=0.9;  %radiator  panel  emissivity 
trad=0.018;  %radiator  panel  thickness  (m) 
wrad=1.5;  %radiator  panel  width  (m) 

Lradl2=9*sin ( 10*pi/180 ) ;  %radiator  panels  1&2  length  (m) 

Aradl2=4*wrad*Lradl2 ;  %radiator  panels  1&2  area  (mA2) 

Lrad34= ( 9+1 . 6 ) *sin ( 10*pi/180) ; %radiator  panels  3&4  length  (m) 

Arad34=4*wrad*Lrad34 ; %radiator  panels  3&4  area  (mA2) 

Lrad56= ( 9+2*1 . 6 ) *sin ( 10*pi/180 ) ; %radiator  panels  5&6  length  (m) 

Arad56=4*wrad*Lrad56 ; %radiator  panels  5&6  area  (mA2) 

Lrad78= ( 9+3*1 . 6 ) *sin ( 10*pi/180 ) ; %radiator  panels  7&8  length  (m) 

Arad78=4*wrad*Lrad78 ; %radiator  panels  7&8  area  (mA2) 

Lrad910= ( 9+4*1 . 6 ) *sin ( 10*pi/180 ) ; %radiator  panels  9&10  length  (m) 

Arad910=4*wrad*Lrad910 ; %radiator  panels  9&10  area  (mA2) 

Lradlll2= ( 9+5*1 . 6 ) *sin ( 10*pi/180 ) ; %radiator  panels  11&12  length  (m) 
Aradlll2=4*wrad*Lradlll2; %radiator  panels  11&12  area  (mA2) 

Lradl314= ( 9+6*1 . 6 ) *sin ( 10*pi/180 ) ; %radiator  panels  13&14  length  (m) 
Aradl314=4*wrad*Lradl314; %radiator  panels  13&14  area  (mA2) 

Lradl516= ( 9+7*1 . 6 ) *sin ( 10*pi/180 ) ; %radiator  panels  15&16  length  (m) 
Aradl516=4*wrad*Lradl516 ; %radiator  panels  15&16  area  (mA2) 

Lradl718= (9+8*1 . 6) *sin (10*pi/180) ;%radiator  panels  17&18  length  (m) 
Aradl718=4*wrad*Lradl718; %radiator  panels  17&18  area  (mA2) 
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Lradl920= ( 9+9*1 . 6 ) *sin ( 10*pi/180 ) ; %radiator  panels  19&20  length  (m) 
Aradl920=4*wrad*Lradl920 ; %radiator  panels  19&20  area  (mA2) 

Lrad2122= ( 9+10*1 . 6 ) *sin ( 10*pi/180 ); %radiator  panels  21&22  length  (m) 
Arad2122=4*wrad*Lrad2122; %radiator  panels  21&22  area  (mA2) 

Lrad2324= ( 9+11*1 . 6 ) *sin ( 10*pi/180 ); %radiator  panels  23&24  length  (m) 
Arad2324=4*wrad*Lrad2324 ; %radiator  panels  23&24  area  (mA2) 

Lrad2526= ( 9+12*1 . 6 ) *sin ( 10*pi/180 ); %radiator  panels  25&26  length  (m) 
Arad2526=4*wrad*Lrad2526 ; %radiator  panels  25&26  area  (mA2) 

Lrad2728= (9+13*1 . 6) *sin (10*pi/180) ;%radiator  panels  27&28  length  (m) 
Arad2728=4*wrad*Lrad2728 ; %radiator  panels  27&28  area  (mA2) 

Lrad2930= ( 9+14*1 . 6 ) *sin ( 10*pi/180 ); %radiator  panels  29&30  length  (m) 
Arad2930=4*wrad*Lrad2930 ; %radiator  panels  29&30  area  (mA2) 

Lrad3132= ( 9+15*1 . 6 ) *sin ( 10*pi/180 ); %radiator  panels  31&32  length  (m) 
Arad3132=4*wrad*Lrad3132 ; %radiator  panels  31&32  area  (mA2) 


(m) 


%Heat  pipes  for  panels  3&4  condenser  thermal 

%Heat  pipes  for  panels  5&6  condenser  thermal 

%Heat  pipes  for  panels  7&8  condenser  thermal 

%Heat  pipes  for  panels  9&10  condenser 
;  %Heat  pipes  for  panels  11&12  condenser 


Lrad3334= (9+16*1 . 6) *sin (10*pi/180) ;%radiator  panels  33&34  length 
Arad3334=4*wrad*Lrad3334 ; %radiator  panels  33&34  area  (mA2) 

Lrad3536= (9+5 . 64*1 . 6) *sin (10*pi/180) ;%radiator  panels  35&36  length  (m) 
Arad3536=4*wrad*Lrad3536 ; %radiator  panels  35&36  area  (mA2) 

Aradtot=Aradl2+Arad34+Arad56+Arad78+Arad910+Aradlll2+Aradl314+Aradl516+Aradl718+Aradl920+ 
Arad2122+Arad2324+Arad2526+Arad2728+Arad2930+Arad3132+Arad3334+Arad3536 ; %radiator  panels 
total  area  (mA2) 

Rpcl2=log (dradhp/idradhp) / ( 4*pi*Lradl2*kp) ;  %Heat  pipes  for  panels  1&2  condenser  thermal 
resistance  (K/W) 

Rpc34=log (dradhp/idradhp) / ( 4*pi*Lrad34*kp) 
resistance  (K/W) 

Rpc56=log (dradhp/idradhp) / ( 4*pi*Lrad56*kp) 
resistance  (K/W) 

Rpc78=log (dradhp/idradhp) / ( 4*pi*Lrad78*kp) 
resistance  (K/W) 

Rpc910=log (dradhp/idradhp) / ( 4*pi*Lrad910*kp) ; 
thermal  resistance  (K/W) 

Rpclll2=log (dradhp/idradhp) / ( 4*pi*Lradlll2*kp) ; 
thermal  resistance  (K/W) 

Rpcl314=log (dradhp/idradhp) / ( 4*pi*Lradl314*kp) ; 
thermal  resistance  (K/W) 

Rpcl516=log (dradhp/idradhp) / ( 4*pi*Lradl516*kp) ; 
thermal  resistance  (K/W) 

Rpcl718=log (dradhp/idradhp) / ( 4*pi*Lradl718*kp) ; 
thermal  resistance  (K/W) 

Rpcl920=log (dradhp/idradhp) / ( 4*pi*Lradl920*kp) ; 
thermal  resistance  (K/W) 

Rpc2122=log (dradhp/idradhp) / ( 4*pi*Lrad2122*kp) ; 
thermal  resistance  (K/W) 

Rpc2324=log (dradhp/idradhp) / ( 4*pi*Lrad2324*kp) ; 
thermal  resistance  (K/W) 

Rpc2526=log (dradhp/idradhp) / ( 4*pi*Lrad2526*kp) ; 
thermal  resistance  (K/W) 

Rpc2728=log (dradhp/idradhp) / ( 4*pi*Lrad2728*kp) ; 
thermal  resistance  (K/W) 

Rpc2930=log (dradhp/idradhp) / ( 4*pi*Lrad2930*kp) ; 
thermal  resistance  (K/W) 

Rpc3132=log (dradhp/idradhp) / ( 4*pi*Lrad3132*kp) ; 
thermal  resistance  (K/W) 

Rpc3334=log (dradhp/idradhp) / ( 4*pi*Lrad3334*kp) ; 
thermal  resistance  (K/W) 

Rpc3536=log (dradhp/idradhp) / ( 4*pi*Lrad3536*kp) ; 
thermal  resistance  (K/W) 

Rpwcl2=log ( idradhp/idhpw) / ( 4*pi*Lradl2*kef f rhp) 
thermal  resistance  (K/W) 

Rpwc34=log ( idradhp/idhpw) / ( 4*pi*Lrad34*kef f rhp) 
thermal  resistance  (K/W) 

Rpwc56=log ( idradhp/idhpw) / ( 4*pi*Lrad56*kef f rhp) 
thermal  resistance  (K/W) 

Rpwc78=log ( idradhp/idhpw) / ( 4*pi*Lrad78*kef f rhp) 


%Heat  pipes  for  panel? 

3  13&14  condenser 

%Heat 

pipes 

for 

panels 

15&16 

condenser 

%Heat 

pipes 

for 

panels 

17&18 

condenser 

%Heat 

pipes 

for 

panels 

19&20 

condenser 

%Heat 

pipes 

for 

panels 

21&22 

condenser 

%Heat 

pipes 

for 

panels 

23&24 

condenser 

%Heat 

pipes 

for 

panels 

25&26 

condenser 

%Heat 

pipes 

for 

panels 

27&28 

condenser 

%Heat 

pipes 

for 

panels 

29&30 

condenser 

%Heat 

pipes 

for 

panels 

31&32 

condenser 

%Heat 

pipes 

for 

panels 

33&34 

condenser 

%Heat 

pipes 

for 

panels 

35&36 

condenser 

;%Heat  pipes  for  panels  1&2  condenser  wick 
;%Heat  pipes  for  panels  3&4  condenser  wick 
; %Heat  pipes  for  panels  5&6  condenser  wick 
;%Heat  pipes  for  panels  7&8  condenser  wick 


thermal  resistance  (K/W) 

Rpwc910=log ( idradhp/idhpw) /( 4*pi*Lrad910*keff rhp) ; %Heat  pipes  for  panels  9&10  condenser 
wick  thermal  resistance  (K/W) 

Rpwclll2=log ( idradhp/idhpw) /( 4*pi*Lradlll2*keff rhp) ; %Heat  pipes  for  panels  11&12 
condenser  wick  thermal  resistance  (K/W) 

Rpwcl314=log ( idradhp/idhpw) /( 4*pi*Lradl314*kef f rhp) ; %Heat  pipes  for  panels  13&14 
condenser  wick  thermal  resistance  (K/W) 
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Rpwcl516=log ( idradhp/idhpw) / ( 4*pi*Lradl516*kef f rhp) ; %Heat  pipes  for  panels  15&16 
condenser  wick  thermal  resistance  (K/W) 

Rpwcl718=log (idradhp/idhpw) / ( 4*pi*Lradl718*kef frhp) ; %Heat  pipes  for  panels  17&18 
condenser  wick  thermal  resistance  (K/W) 

Rpwcl920=log ( idradhp/idhpw) /( 4*pi*Lradl920*kef frhp) ; %Heat  pipes  for  panels  19&20 
condenser  wick  thermal  resistance  (K/W) 

Rpwc2122=log ( idradhp/idhpw) /( 4*pi*Lrad2122*kef frhp) ; %Heat  pipes  for  panels  21&22 
condenser  wick  thermal  resistance  (K/W) 

Rpwc2324=log ( idradhp/idhpw) /( 4*pi*Lrad2324*kef frhp) ; %Heat  pipes  for  panels  23&24 
condenser  wick  thermal  resistance  (K/W) 

Rpwc2526=log ( idradhp/idhpw) /( 4*pi*Lrad2526*kef frhp) ; %Heat  pipes  for  panels  25&26 
condenser  wick  thermal  resistance  (K/W) 

Rpwc2728=log ( idradhp/idhpw) /( 4*pi*Lrad2728*kef frhp) ; %Heat  pipes  for  panels  27&28 
condenser  wick  thermal  resistance  (K/W) 

Rpwc2930=log ( idradhp/idhpw) /( 4*pi*Lrad2930*kef frhp) ; %Heat  pipes  for  panels  29&30 
condenser  wick  thermal  resistance  (K/W) 

Rpwc3132=log ( idradhp/idhpw) /( 4*pi*Lrad3132*kef frhp) ; %Heat  pipes  for  panels  31&32 
condenser  wick  thermal  resistance  (K/W) 

Rpwc3334=log ( idradhp/idhpw) /( 4*pi*Lrad3334*kef frhp) ; %Heat  pipes  for  panels  33&34 
condenser  wick  thermal  resistance  (K/W) 

Rpwc3536=log ( idradhp/idhpw) /( 4*pi*Lrad3536*kef frhp) ; %Heat  pipes  for  panels  35&36 
condenser  wick  thermal  resistance  (K/W) 

Rpel2=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ; %Heat  pipes  for  panels  1&2  evaporator  thermal 
resistance  (K/W) 

Rpe34=log (dradhp/idradhp) /( 4*pi*Lerhp*kp) ; %Heat  pipes  for  panels  3&4  evaporator  thermal 
resistance  (K/W) 

Rpe56=log (dradhp/idradhp) /( 4*pi*Lerhp*kp) ; %Heat  pipes  for  panels  5&6  evaporator  thermal 
resistance  (K/W) 

Rpe78=log (dradhp/idradhp) /( 4*pi*Lerhp*kp) ; %Heat  pipes  for  panels  7&8  evaporator  thermal 
resistance  (K/W) 

Rpe910=log (dradhp/idradhp) /( 4*pi*Lerhp*kp) ; %Heat  pipes  for  panels  9&10  evaporator  thermal 
resistance  (K/W) 


Rpelll2=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

11&12 

evaporator 

Rpel314=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

13&14 

evaporator 

Rpel516=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

15&16 

evaporator 

Rpel718=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

17&18 

evaporator 

Rpel920=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

19&20 

evaporator 

Rpe2122=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

21&22 

evaporator 

Rpe2324=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

23&24 

evaporator 

Rpe2526=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

25&26 

evaporator 

Rpe2728=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

27&28 

evaporator 

Rpe2930=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

29&30 

evaporator 

Rpe3132=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

31&32 

evaporator 

Rpe3334=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

33&34 

evaporator 

Rpe3536=log (dradhp/idradhp) / ( 4*pi*Lerhp*kp) ;%Heat 
thermal  resistance  (K/W) 

pipes 

for 

panels 

35&36 

evaporator 

Rpwel2=log ( idradhp/idhpw) /( 4*pi*Lerhp*kef frhp) ; %Heat  pipes  for  panels  1&2  evaporator  wick 
thermal  resistance  (K/W) 

Rpwe34=log ( idradhp/idhpw) /( 4*pi*Lerhp*kef frhp) ; %Heat  pipes  for  panels  3&4  evaporator  wick 
thermal  resistance  (K/W) 

Rpwe56=log ( idradhp/idhpw) /( 4*pi*Lerhp*kef frhp) ; %Heat  pipes  for  panels  5&6  evaporator  wick 
thermal  resistance  (K/W) 

Rpwe78=log ( idradhp/idhpw) /( 4*pi*Lerhp*kef frhp) ; %Heat  pipes  for  panels  7&8  evaporator  wick 
thermal  resistance  (K/W) 

Rpwe910=log ( idradhp/idhpw) /( 4*pi*Lerhp*kef frhp) ; %Heat  pipes  for  panels  9&10  evaporator 
wick  thermal  resistance  (K/W) 

Rpwelll2=log ( idradhp/idhpw) /( 4*pi*Lerhp*kef frhp) ; %Heat  pipes  for  panels  11&12  evaporator 
wick  thermal  resistance  (K/W) 
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%Heat 

pipes 

for 

panels 

13&14 

evaporator 

%Heat 

pipes 

for 

panels 

15&16 

evaporator 

%Heat 

pipes 

for 

panels 

17&18 

evaporator 

%Heat 

pipes 

for 

panels 

19&20 

evaporator 

%Heat 

pipes 

for 

panels 

21&22 

evaporator 

%Heat 

pipes 

for 

panels 

23&24 

evaporator 

%Heat 

pipes 

for 

panels 

25&26 

evaporator 

%Heat 

pipes 

for 

panels 

27&28 

evaporator 

%Heat 

pipes 

for 

panels 

29&30 

evaporator 

%Heat 

pipes 

for 

panels 

31&32 

evaporator 

%Heat 

pipes 

for 

panels 

33&34 

evaporator 

%Heat 

pipes 

for 

panels 

35&36 

evaporator 

wick  thermal  resistance  (K/W) 

Rpwel516=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rpwel718=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rpwel920=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rpwe2122=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rpwe2324=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rpwe2526=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rpwe2728=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rpwe2930=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rpwe3132=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rpwe3334=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rpwe3536=log ( idradhp/idhpw) / ( 4*pi*Lerhp*kef f rl 
wick  thermal  resistance  (K/W) 

Rtotl2=l/ ( 2*hradhp*Aerhp) +Rpel2+Rpwel2+Rpcl2+Rpwcl2 ; %Heat  pipes  for  panels  1&2  total 
thermal  resistance  (K/W) 

Rtot34=l/ ( 2*hradhp*Aerhp) +Rpe34+Rpwe34+Rpc34+Rpwc34; %Heat  pipes  for  panels  3&4  total 
thermal  resistance  (K/W) 

Rtot56=l/ ( 2*hradhp*Aerhp) +Rpe56+Rpwe56+Rpc56+Rpwc56 ; %Heat  pipes  for  panels  5&6  total 
thermal  resistance  (K/W) 

Rtot78=l/ ( 2*hradhp*Aerhp) +Rpe78+Rpwe78+Rpc78+Rpwc78 ; %Heat  pipes  for  panels  7&8  total 
thermal  resistance  (K/W) 

Rtot910=l/ ( 2*hradhp*Aerhp) +Rpe910+Rpwe910+Rpc910+Rpwc910 ; %Heat  pipes  for  panels  9&10 
total  thermal  resistance  (K/W) 

Rtotlll2=l/ ( 2*hradhp*Aerhp) +Rpelll2+Rpwelll2+Rpclll2+Rpwclll2 ; %Heat  pipes  for  panels 
11&12  total  thermal  resistance  (K/W) 

Rtotl314=l/ ( 2*hradhp*Aerhp) +Rpel314+Rpwel314+Rpcl314+Rpwcl314 ; %Heat  pipes  for  panels 
13&14  total  thermal  resistance  (K/W) 

Rtotl516=l/ ( 2*hradhp*Aerhp) +Rpel516+Rpwel516+Rpcl516+Rpwcl516 ; %Heat  pipes  for  panels 
15&16  total  thermal  resistance  (K/W) 

Rtotl718=l/ ( 2*hradhp*Aerhp) +Rpel718+Rpwel718+Rpcl718+Rpwcl718 ; %Heat  pipes  for  panels 
17&18  total  thermal  resistance  (K/W) 

Rtotl920=l/ ( 2*hradhp*Aerhp) +Rpel920+Rpwel920+Rpcl920+Rpwcl920 ; %Heat  pipes  for  panels 
19&20  total  thermal  resistance  (K/W) 

Rtot2122=l/ ( 2*hradhp*Aerhp) +Rpe2122+Rpwe2122+Rpc2122+Rpwc2122 ; %Heat  pipes  for  panels 
21&22  total  thermal  resistance  (K/W) 

Rtot2324=l/ ( 2*hradhp*Aerhp) +Rpe2324+Rpwe2324+Rpc2324+Rpwc2324 ; %Heat  pipes  for  panels 
23&24  total  thermal  resistance  (K/W) 

Rtot2526=l/ ( 2*hradhp*Aerhp) +Rpe2526+Rpwe2526+Rpc2526+Rpwc2526 ; %Heat  pipes  for  panels 
25&26  total  thermal  resistance  (K/W) 

Rtot2728=l/ ( 2*hradhp*Aerhp) +Rpe2728+Rpwe2728+Rpc2728+Rpwc2728 ; %Heat  pipes  for  panels 
27&28  total  thermal  resistance  (K/W) 

Rtot2930=l/ ( 2*hradhp*Aerhp) +Rpe2930+Rpwe2930+Rpc2930+Rpwc2930 ; %Heat  pipes  for  panels 
29&30  total  thermal  resistance  (K/W) 

Rtot3132=l/ ( 2*hradhp*Aerhp) +Rpe3132+Rpwe3132+Rpc3132+Rpwc3132 ; %Heat  pipes  for  panels 
31&32  total  thermal  resistance  (K/W) 

Rtot3334=l/ ( 2*hradhp*Aerhp) +Rpe3334+Rpwe3334+Rpc3334+Rpwc3334 ; %Heat  pipes  for  panels 
33&34  total  thermal  resistance  (K/W) 

Rtot3536=l/ ( 2*hradhp*Aerhp) +Rpe3536+Rpwe3536+Rpc3536+Rpwc3536 ; %Heat  pipes  for  panels 
35&36  total  thermal  resistance  (K/W) 
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17&18 

(K) 

Initial 

flow 

temperature 

before 

heat 

pipes 

19&20 

(K) 

Initial 

flow 

temperature 

before 

heat 

pipes 

21&22 

(K) 

Initial 

flow 

temperature 

before 

heat 

pipes 

23&24 

(K) 

Initial 

flow 

temperature 

before 

heat 

pipes 

25&26 

(K) 

Initial 

flow 

temperature 

before 

heat 

pipes 

27&28 

(K) 

Initial 

flow 

temperature 

before 

heat 

pipes 

29&30 

(K) 

Initial 

flow 

temperature 

before 

heat 

pipes 

31&32 

(K) 

Initial 

flow 

temperature 

before 

heat 

pipes 

33&34 

(K) 

Initial 

flow 

temperature 

before 

heat 

pipes 

35&36 

(K) 

Initial 

flow 

temperature 

after  heat  pipes  : 

35&36  i 

(K) 

Tx=427 . 730225305744; %Initial  flow  temperature  before  the  NaK  channels  (K) 
Tabrl= (Tarl+Tbrl) /2; %Initial  average  flow  temperature  at  heat  pipes  1&2  (K) 

Tbcrl= (Tbrl+Tcrl) /2; %Initial  average  flow  temperature  at  heat  pipes  3&4  (K) 

Tcdrl= (Tcrl+Tdrl) /2; %Initial  average  flow  temperature  at  heat  pipes  5&6  (K) 

Tderl= (Tdrl+Terl) /2; %Initial  average  flow  temperature  at  heat  pipes  7&8  (K) 
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(K) 

Initial 
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19&20 

(K) 

Initial 
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flow 

temperature 
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21&22 

(K) 

Initial 

average 

flow 

temperature 

at 

heat 

pipes 

23&24 

(K) 
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average 
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temperature 
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pipes 

25&26 

(K) 
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average 
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temperature 
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heat 

pipes 

27&28 

(K) 

Initial 

average 

flow 

temperature 

at 
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pipes 

29&30 

(K) 

Initial 

average 

flow 

temperature 

at 
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31&32 

(K) 
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flow 

temperature 
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pipes 

33&34 

(K) 

Initial 

average 

flow 

temperature 

at 

heat 

pipes 

35&36 

(K) 

itial  average  flow  temperature  at  NaK  channels 

(K) 

dQpl2e=mdott*cptl* (Tarl-Tbrl) ; %heat  transfer  into  heat  pipes  1&2  (W) 
dQp34e=mdott*cptl* (Tbrl-Tcrl) ; %heat  transfer  into  heat  pipes  3&4  (W) 
dQp56e=mdott*cptl* (Tcrl-Tdrl) ; %heat  transfer  into  heat  pipes  5&6  (W) 
dQp78e=mdott*cptl* (Tdrl-Terl) ; %heat  transfer  into  heat  pipes  7&8  (W) 


dQplll2e=mdott*cptl* (Tfrl-Tgrl) 
dQpl314e=mdott*cptl* (Tgrl-Thrl) 
dQpl516e=mdott*cptl* (Thrl-Tirl) 
dQpl718e=mdott*cptl* (Tirl-Tjrl) 
dQpl920e=mdott*cptl* (Tjrl-Tkrl) 
dQp2122e=mdott*cptl* (Tkrl-Tlrl) 
dQp2324e=mdott*cptl* (Tlrl-Tmrl) 
dQp2526e=mdott*cptl* (Tmrl-Tnrl) 
dQp2728e=mdott*cptl* (Tnrl-Torl) 
dQp2930e=mdott*cptl* (Torl-Tprl) 
dQp3132e=mdott*cptl* (Tprl-Tqrl) 
dQp3334e=mdott*cptl* (Tqrl-Trrl) 
dQp3536e=mdott*cptl* (Trrl-Tsrl) 

Tabp=Tabrl-dQpl2e*Rtotl2 ; %temperature  of  radiator  panels 
Tbcp=Tbcrl-dQp34e*Rtot34; %temperature  of  radiator  panels 
Tcdp=Tcdrl-dQp56e*Rtot56 ; %temperature  of  radiator  panels 
Tdep=Tderl-dQp78e*Rtot78 ; %temperature  of  radiator  panels 
Tefp=Tef-dQp910e*Rtot910 ; %temperature  of  radiator  panels 


;heat  transfer  into  heat  pipes  9&10  (W) 

%heat 

transfer 

into 

heat 

pipes 

11&12 

(W) 

%heat 

transfer 

into 

heat 

pipes 

13&14 

(W) 

%heat 

transfer 

into 

heat 

pipes 

15&16 

(W) 

%heat 

transfer 

into 

heat 

pipes 

17&18 

(W) 

%heat 

transfer 

into 

heat 

pipes 

19&20 

(W) 

%heat 

transfer 

into 

heat 

pipes 

21&22 

(W) 

%heat 

transfer 

into 

heat 

pipes 

23&24 

(W) 

%heat 

transfer 

into 

heat 

pipes 

25&26 

(W) 

%heat 

transfer 

into 

heat 

pipes 

27&28 

(W) 

%heat 

transfer 

into 

heat 

pipes 

29&30 

(W) 

%heat 

transfer 

into 

heat 

pipes 

31&32 

(W) 

%heat 

transfer 

into 

heat 

pipes 

33&34 

(W) 

%heat 

transfer 

into 

heat 

pipes 

35&36 

(W) 

1&2 

3&4 

5&6 

7&8 


(K) 

(K) 

(K) 

(K) 


9&10  (K) 


Ti jp=Ti j-dQpl718e* 


%temperature 

of 

radiator 

panels 

11&12 

(K) 

%temperature 

of 

radiator 

panels 

13&14 

(K) 

%temperature 

of 

radiator 

panels 

15&16 

(K) 

%temperature 

of 

radiator 

panels 

17&18 

(K) 

%temperature 

of 

radiator 

panels 

19&20 

(K) 

%temperature 

of 

radiator 

panels 

21&22 

(K) 

%temperature 

of 

radiator 

panels 

23&24 

(K) 

%temperature 

of 

radiator 

panels 

25&26 

(K) 

%temperature 

of 

radiator 

panels 

27&28 

(K) 

% temperature 

of 

radiator 

panels 

29&30 

(K) 

%temperature 

of 

radiator 

panels 

31&32 

(K) 

%temperature 

of 

radiator 

panels 

33&34 

(K) 

%temperature 
erature  (K) 

of 

radiator 

panels 

35&36 

(K) 

(TabpA4-Tamb/ 

'4 )  ; 

:  %heat  rejected 

by  panels  1&2  (W) 

Tamb=4;  %ambient  space  t 
dQpl2=emiss*sigmasb*Arad 

dQp34=emiss*sigmasb*Arad34*  (Tbcp/'4-TambA4 )  , 
dQp56=emiss*sigmasb*Arad56*  (Tcdp/'4-TambA4 )  , 
dQp78=emiss*sigmasb*Arad78*  (Tdep/'4-TambA4 )  , 
dQp910=emiss*sigmasb*Arad910*  (Tefp/'4-TambA4)  ; %heat  rejected  by  panels  9&10  (W) 
dQplll2=emiss*sigmasb*Aradlll2* (TfgpA4-TambA4) ; %heat  rejected  by  panels  11&12  (W) 


;%heat  rejected  by  panels  3&4  (W) 
;%heat  rejected  by  panels  5&6  (W) 
;%heat  rejected  by  panels  7&8  (W) 


F10 


iheat 

rejected 

by 

panels 

13&14 

(W) 

iheat 

rejected 

by 

panels 

15&16 

(W) 

iheat 

rejected 

by 

panels 

17&18 

(W) 

iheat 

rejected 

by 

panels 

19&20 

(W) 

iheat 

rejected 

by 

panels 

21&22 

(W) 

iheat 

rejected 

by 

panels 

23&24 

(W) 

iheat 

rejected 

by 

panels 

25&26 

(W) 

iheat 

rejected 

by 

panels 

27&28 

(W) 

iheat 

rejected 

by 

panels 

29&30 

(W) 

iheat 

rejected 

by 

panels 

31&32 

(W) 

iheat 

rejected 

by 

panels 

33&34 

(W) 

!*pi) 1 

k (TrspA4-TambA4) )  ;  %. 

length 

of 

3&4 

5&6 

7&8 


(mA3 ) 
(mA3 ) 
(mA3 ) 
(mA3 ) 


dQpl314=emiss*sigmasb*Aradl314* (TghpA4-TambA' 
dQpl516=emiss*sigmasb*Aradl516* (ThipA4-TambA' 
dQpl718=emiss*sigmasb*Aradl718* (Ti jpA4-TambA' 
dQpl920=emiss*sigmasb*Aradl920* (T jkpA4-TambA 
dQp2122=emiss*sigmasb*Arad2122* (TklpA4-TambA‘ 
dQp2324=emiss*sigmasb*Arad2324* (TlmpA4-TambA‘ 
dQp2526=emiss*sigmasb*Arad2526* (TmnpA4-TambA‘ 
dQp2728=emiss*sigmasb*Arad2728* (TnopA4-TambA‘ 
dQp2930=emiss*sigmasb*Arad2930* (ToppA4-TambA 
dQp3132=emiss*sigmasb*Arad3132* (TpqpA4-TambA‘ 
dQp3334=emiss*sigmasb*Arad3334* (TqrpA4-TambA' 

Lrad3536=dQp3536e/ ( 4*emiss*sigmasb*wrad*sin ( : 
radiators  35&36  (m) 

Arad3536=4*wrad*Lrad3536*sin ( 10*pi/180) ; %area  of  radiators  35&36  (mA2) 
dQp3536=emiss*sigmasb*Arad3536* (TrspA4-TambA4) ; %heat  rejection  of  radiators  35&36  (W) 
dQpt  o t =dQp 1 2  +dQp3  4 +dQp  5  6 +dQp  7  8 +dQp  910 +dQp 1112 +dQp 1314 +dQp 1516  +dQp 1718  +dQp 1920  +dQp2 122  +dQp 
2324+dQp2526+dQp2728+dQp2930+dQp3132+dQp3334+dQp3536 ; %total  radiator  heat  rejection  (W) 
wsad=0.0032;  %graphite  heat  pipe  saddle  thickness  (m) 
ssad=dradhp+2*wsad;  %graphite  heat  pipe  saddle  outer  diameter  (m) 
psad=550;  %graphite  heat  pipe  density  (kg/mA3) 
pfin=180;  %radiator  panel  density  (kg/mA3) 

cpg=720;  %graphite  heat  pipe  saddle  specific  heat  ( J/ (kg*K) ) 

Vradl2=Aradl2*trad+Lradl2* ( ssadA2-pi*dradhpA2/4 ) ; %volume  of  radiators  1&2 
Vrad34=Arad34*trad+Lrad34* (ssadA2-pi*dradhpA2/4) ;%volume  of  radiators 
Vrad56=Arad56*trad+Lrad56* ( ssadA2-pi*dradhpA2/4 ) ;%volume  of  radiators 
Vrad78=Arad78*trad+Lrad78* (ssadA2-pi*dradhpA2/4) ;%volume  of  radiators 
Vrad910=Arad910*trad+Lrad910* (ssadA2-pi*dradhpA2/4) ; %volume  of  radiators  9&10  (mA3) 
Vradlll2=Aradlll2*trad+Lradlll2* ( ssadA2-pi*dradhpA2/4 ) ; %volume  of  radiators  11&12  (mA3) 
Vradl314=Aradl314*trad+Lradl314* ( ssadA2-pi*dradhpA2/4 ) ; %volume  of  radiators  13&14 
Vradl516=Aradl516*trad+Lradl516* ( ssadA2-pi*dradhpA2/4 ) ; %volume  of  radiators  15&16 
Vradl718=Aradl718*trad+Lradl718* (ssadA2-pi*dradhpA2/4) ; %volume  of  radiators  17&18 
Vradl920=Aradl920*trad+Lradl920* ( ssadA2-pi*dradhpA2/4 ) ; %volume  of  radiators  19&20 
Vrad2122=Arad2122*trad+Lrad2122* ( ssadA2-pi*dradhpA2/4 ) ; %volume  of  radiators  21&22 
Vrad2324=Arad2324*trad+Lrad2324* ( ssadA2-pi*dradhpA2/4 ) ; %volume  of  radiators  23&24 
Vrad2526=Arad2526*trad+Lrad2526* (ssadA2— pi*dradhpA2/4) ; %volume  of  radiators  25&26 
Vrad2728=Arad2728*trad+Lrad2728* (ssadA2-pi*dradhpA2/4) ;%volume  of  radiators 
Vrad2930=Arad2930*trad+Lrad2930* ( ssadA2-pi*dradhpA2/4 ) ;%volume  of  radiators 
Vrad3132=Arad3132*trad+Lrad3132* ( ssadA2-pi*dradhpA2/4 ) ; %volume  of  radiators  31&32 
Vrad3334=Arad3334*trad+Lrad3334* ( ssadA2-pi*dradhpA2/4 ) ; %volume  of  radiators  33&34 
Vrad3536=Arad3536*trad+Lrad3536* (ssadA2-pi*dradhpA2/4) ; %volume  of  radiators  35&36 
mradl2=pf in*Aradl2*trad+psad*Lradl2* ( ssadA2-pi*dradhpA2/4 ) ; %mass  of  radiators  1&2 
mrad34=pf in*Arad34*trad+psad*Lrad34* ( ssadA2-pi*dradhpA2/4 ) ; %mass  of  radiators  3&4 
mrad56=pf in*Arad56*trad+psad*Lrad56* ( ssadA2-pi*dradhpA2/4 ) ; %mass  of  radiators  5&6 
mrad78=pf in*Arad78*trad+psad*Lrad78* ( ssadA2-pi*dradhpA2/4 ) ; %mass  of  radiators  7&8 
mrad910=pf in*Arad910*trad+psad*Lrad910* ( ssadA2-pi*dradhpA2/4 ) ; %mass  of  radiators  9&10 
(mA3 ) 

mradlll2=pf in*Aradlll2*trad+psad*Lradlll2* ( ssadA2-pi*dradhpA2/4 ) ; %mass  of  radiators  11&12 
(mA3 ) 

mradl314=pf in*Aradl314*trad+psad*Lradl314* (ssadA2-pi*dradhpA2/4) ; %mass  of  radiators  13&14 
(mA3 ) 

mradl516=pf in*Aradl516*trad+psad*Lradl516* (ssadA2-pi*dradhpA2/4) ; %mass  of  radiators  15&16 
(mA3 ) 

mradl718=pf in*Aradl718*trad+psad*Lradl718* (ssadA2-pi*dradhpA2/4) ; %mass  of  radiators  17&18 
(mA3 ) 

mradl920=pf in*Aradl920*trad+psad*Lradl920* ( ssadA2-pi*dradhpA2/4 ) ;%mass  of  radiators 
(mA3 ) 

mrad2122=pf in*Arad2122*trad+psad*Lrad2122* ( ssadA2-pi*dradhpA2/4 ) ; %mass  of  radiators  21&22 
(mA3 ) 

mrad2324=pf in*Arad2324*trad+psad*Lrad2324* (ssadA2- 
mrad2526=pf in*Arad2526  *trad+psad*Lrad2526  * (ssadA2- 
mrad2728=pf in*Arad2728*trad+psad*Lrad2728* (ssadA2- 
mrad2930=pf in*Arad2930*trad+psad*Lrad2930* ( ssadA2- 
(mA3 ) 

mrad3132=pf in*Arad3132*trad+psad*Lrad3132* (ssadA2-pi*dradhpA2/4) ; %mass  of  radiators  31&32 
(mA3 ) 

mrad3334=pf in*Arad3334*trad+psad*Lrad3334* ( ssadA2-pi*dradhpA2/4 ) ;%mass  of  radiators 
(mA3 ) 

mrad3536=pf in*Arad3536*trad+psad*Lrad3536* (ssadA2-pi*dradhpA2/4) ;%mass  of  radiators 
(mA3 ) 

Vhpx=Arlp*dradhp-2*drlp*pi*dradhpA2/4 ;  %flow  volume  at  radiator  heat  pipe  condensers 
(mA3) 


27&28 

29&30 


(mA3 ) 
(mA3) 
(mA3) 
(mA3) 
(mA3) 
(mA3) 
(mA3) 
(mA3 ) 
(mA3) 
(mA3) 
(mA3) 
(mA3) 
(mA3) 
(mA3) 
(mA3 ) 
(mA3) 
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of 
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of 
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13&14 

of 
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of 
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of 

radiators 

19&20 

of 

radiators 

21&22 

of 

radiators 

23&24 

of 

radiators 

25&26 

of 

radiators 

27&28 

of 

radiators 

29&30 

of 

radiators 

31&32 

of 

radiators 

33&34 

of 
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35&36 

pe 
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n 
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To  Workspace4 

X 
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Ff 
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Embedded  MATLAB  Function 
function 

dTHx=LMTD ( Tbl , Tb2 , Tc2 , Tel ) 
if  abs ( (Tbl-Tc2 ) - (Tb2-Tcl ) ) <=1 
dTHx=Tbl-Tc2 ; 

else 

dTHx= (Tbl-Tc2- (Tb2- 
Tel ) ) /log ( (Tbl-Tc2) / (Tb2-Tcl) ) ; 
end 


Stirling  Loop 
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Freq 

function  [freq,t2]  = 

fen (dV, t, tl, freql ) 

f req=f reql ; 

t2=tl ; 

if  abs (dV) <=1*10A (-2) 
freq=l/ (t-t2 ) ; 
t2=t ; 

end 

if  abs (t-tl) <=1/ (2*95) 
f req=f reql ; 
t2=tl ; 

end 

if  freq  >=  130 
f req=f req/ 2 ; 

end 


dQk 

function  [dQk,Qk,t2]  = 

fen (dQkO, dQkl, t, tl, dQk 2, Qkl, Qk2) 

dQk=dQk2 ; 

Qk=Qkl; 

t2=tl; 

if  dQk0>0  &&  dQklcO 

dQk= (Qk2-Qkl ) / (t-t2) ; 

Qk=Qk2 ; 
t2=t  ; 

end 


dQh 

function  [dQh,Qh,t2]  = 

fen (dQhO , dQhl , t , tl , dQh2 , Qhl , Qh2 ) 

dQh=dQh2 ; 

Qh=Qhl ; 
t2=tl; 

if  dQh0>0  &&  dQhl<0 

dQh= (Qh2-Qhl ) / (t-t2) ; 

Qh=Qh2 ; 
t2=t  ; 

end 
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G10 


Stirling  Main  Block 

function  [pcs, pe, dQk, dQr , dQh, dmcs, dmk, dmr , dmh, dme, hAk, hAr , hAh, Th, Tk]  = 

fen (dmeso, peo, mwork, dVcs, dVe,  Vcs,  Ve,  Twho,  Tmro,  Twko,  t,  dQho,  dQro,  dQko,  hAko,  hAro 

, hAho ) 

dmcs=dmcso;  %initial  compression  space  mass  flow  rate  (kg/s) 

pe=peo;  %initial  expansion  space  mass  pressure  (Pa) 

vc=l . 3 ;  %compression  volume  factor 

ve=l;  %expansion  space  volume  factor 

freq=90;  %stirling  frequency  (Hz) 

Ap=0.0381;  %piston  area  (mA2) 
xpo=0;  %initial  piston  position  (m) 
vpo=0 . 013*2*pi*freq; 

kalt=38.5;  %alternator  constant  (Vs/m) 

Ralt=0.028;  %alternator  resistance  (ohm) 

Halt=0.0002;  %alternator  inductance  (H) 

Zalt=sqrt (RaltA2+ (2*pi*freq*Halt) A2 ) ;  %alternator  impedance  (ohms) 
Ialt=kalt*vpo/Zalt;  %alternator  current  (A) 

Rload=0 . 6835;  %alternator  load  resistance  (ohms) 

Dp=kaltA2/Zalt;  %piston  damping  (kg/s) 
mp=25;  %piston  mass  (kg) 
md=8.5;  %displacer  mass  (kg) 
drd=0.032;  %displacer  rod  diameter  (m) 

Arod=pi*drdA2/4 ;  %displacer  rod  area  (mA2) 

Ade=0. 03892;  %displacer  area-exp.  side  (mA2) 

Adc=Ap;  %displacer  area-comp,  side  (mA2) 
phase=60;  %piston-displacer  phase  angle  (deg) 

xdo=0 . 0084*sin (pi*phase/180 ) ; %displacer  initial  position  (m) 

vdo=0 . 0082*2*pi*freq*cos (pi*phase/180 ) ;  %initial  displacer  velocity  (m/s) 

Dd=0;  %displacer  damping  (kg/s) 

kpl=0;  %lst  order  piston  spring  constant  (N/m) 

kd3=0;  %3rd  order  piston  spring  constant  (N/m) 

kd=3.1*10A6;  %displacer  spring  constant  (Pa) 

pmean=13 . 6*10 A6 ;  %piston  mean  pressure  (Pa) 

Twk=Twko;  %cooler  temperature  (K) 

Twh=Twho;  %heater  temperature  (K) 

Twr= (Twh+Twk) /2;  %regen .  temperature  (K) 

Tmr=Tmro;  %regen.  matrix  temp.  (K) 

Vcso=vc*6*10A (-4)  ;  %compression  space  volume  (mA3) 

Nkchan=3100*4;  %number  cold  exchanger  channels  (dl) 

CHxLg=0.075;  %cold  exch.  length  (m) 
dtk=0 . 0014/2 ;  %cooler  channel  diameter  (m) 

CHxpw= (pi*dtk ) *Nkchan;  %cooler  wetted  perimeter  (m) 

Awk=CHxpw*CHxLg;  %cooler  wall  area  (mA2) 

ACHx=Nkchan*pi*dtkA2/4;  %cooler  cross-sectional  area  (mA2) 

CHxDh=4*ACHx/CHxpw;  %cooler  hydraulic  diameter  (m) 

Vk=ACHx*CHxLg;  %cooler  volume  (mA3) 
rLg=0.055;  %regen.  length  (m) 
por=0.88;  %regen.  porosity  (dl) 

Acr=5 . 017*10A (-2 ) ;  %regen.  area  (mA2) 

Ar=Acr*por;  %regen.  Flow  area  (mA2) 
dr=sqrt ( 4*Ar /pi ) ;  %regen.  diameter  (m) 

Awr=pi*dr*rLg;  %regen.  wall  area  (mA2) 
rwd=25*10 A (-6 ) ;  %regen.  wire  diameter  (m) 
rpw=4* (1-por) *Acr/rwd;  %reg.  wetted  perimeter  (m) 
rDh=4*Ar /rpw;  %regen.  hyd .  diameter  (m) 

Amr=rpw*rLg;  %regen.  matrix  area  (mA2) 
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Gil 


Vr=rLg*Ar;  %regen.  volume  (mA3) 

HHxLg=0.075;  %heater  length  (m) 
htd=. 0013/2;  %heater  tube  diameter  (m) 

Nht=2750*4;  %number  heater  tubes  (dl) 

AHHx=Nht*pi/4*htdA2 ;  %heater  cross-sectional  area  (mA2) 
Vh=AHHx*HHxLg;  %heater  volume  (mA3) 

HHxpw=Nht*pi*htd;  %heater  wetted  perimeter  (m) 

HHxDh=4*AHHx/HHxpw;  %heater  diameter  hyd .  (m) 

Awh=HHxpw*HHxLg;  %heater  wall  area  (mA2) 

Veo=ve*4 . 34*10A (-4 ) ;  %expansion  space  volume  (mA3) 

Vbso=14 . 5*10A  (-3 )  ;  %bounce  space  volume  (mA3) 

Vt=Vcso+Veo+Vh+Vr+Vk;  %total  volume  (mA3) 
gammaHe=l . 667;  %helium  ratio  of  specific  heats  (dl) 

Rgas=2077;  %helium  gas  constant  (J/ (kg*K) ) 
kHe=0.14;  %helium  conductivity  (J/ (kg*K) ) 
cp=5230;  %helium  specific  heat-const,  press.  ( J/ (kg*K) ) 
cv=cp-Rgas;  %helium  specific  heat-const,  vol .  (J/ (kg*K) ) 

cs=500;  %steel  specific  heat  (J/ (kg*K) ) 

CSuther land=l . 286e-6 ;  %Suther land ' s  law  constant 
SSutherland=8 . 63 ;  %Suther land ' s  law  constant 
Diss=0;  %initial  total  dissipation  (W) 

Dissk=0;  %initial  cooler  dissipation  (W) 

Dissr=0;  %initial  regen.  dissipation  (W) 

Dissh=0;  %initial  heater  dissipation  (W) 

pcso=pmean;  %initial  compression  space  pressure  (Pa) 
peo=pmean;  %initial  expansion  space  pressure  (Pa) 
pko=pmean;  %initial  cooler  space  pressure  (Pa) 
pro=pmean;  %initial  regen.  space  pressure  (Pa) 
pho=pmean;  %initial  heater  space  pressure  (Pa) 

WpV=53 . 6*10 A3 ;  %expected  cycle  PV  power  (W) 
eff=0.316;  %expected  cycle  efficiency 
dQhss=WpV/ef f ;  %expected  cycle  heat  addition  (W) 
dQkss=WpV-dQhss ;  %expected  cycle  heat  rejection  (W) 
if  t<=0 . 02 

Th=Twh-dQhss/hAho;  %heater  space  temp.  (K) 

Tk=Twk-dQkss/hAko;  %cooler  space  temp.  (K) 

else 

Th=Twh-dQho/hAho ;  %heater  space  temp.  (K) 

Tk=Twk-dQko/hAko;  %cooler  space  temp.  (K) 

end 

Tr= (Th-Tk) /log (Th/Tk) ;  %regen.  1  temp.  (K) 

Tcs=Tk;  %compression  space  temperature  (K) 

Te=Th;  %expansion  space  temperature  (K) 

muk=CSutherland*TkAl . 5/ (Tk+SSutherland) ;  %cooler  viscosity 
muh=CSutherland*ThAl . 5/ (Th+SSutherland)  ;  %heater  viscosity 
mur=CSutherland*Tr A1 . 5/ (Tr+SSutherland) ;  %regen.  viscosity 
meo=peo*Veo/ (Rgas*Th) ;  %initial  expansion  space  mass  (kg) 
mho=pho*Vh/ (Rgas*Th) ;  %initial  heater  space  mass  (kg) 
mro=pro*Vr / (Rgas*Tr ) ;  %initial  regen.  space  mass  (kg) 
mko=pko*Vk/ (Rgas*Tk ) ;  %initial  cooler  space  mass  (kg) 
mcso=pcso*Vcso/ (Rgas*Tk ) ;  %initial  compression  space  mass  (kg) 
mbo=pmean*Vbso/ (Rgas*Tk )  ;  %initial  bounce  space  mass  (kg) 
mwork= (meo+mho+mro+mko+mcso ) ;  %working  space  mass  (kg) 
kHek= (17 . 201+0 . 06764*Tk-l . 25*10 A (-6 ) *TkA2 ) *418 . 68/10 A5;  %cooler  gas 
conductivity  (W/ (m*K) ) 
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kHer= (17 . 201+0 . 06764*Tr-l . 25*10 A (-6 ) *Tr A2 ) *418 . 68/10 A5;  %regen.  gas 
conductivity  (W/ (m*K) ) 

kHeh= (17 . 201+0 . 06764*Th-l . 25*10 A (-6 ) *ThA2 ) *418 . 68/10 A5;  %heater  gas 
conductivity  (W/ (m*K) ) 

Prk=muk*cp/kHek ;  %cooler  gas  Prandtl  number 

Prr=mur *cp/kHer ;  %regen.  gas  Prandtl  number 

Prh=muh*cp/kHeh;  %heater  gas  Prandtl  number 

Kk=1.5;  %cooler  local  loss  coefficient 
Kr=0;  %regen.  local  loss  coefficient 
Kh=1.5;  %heater  local  loss  coefficient 

pcs=mwork*Rgas/ (Vcs/Tk+Vk/Tk+Vr/Tr+Vh/Th+Ve/Th) ;  %working  space  pressure  (Pa) 
dpcs=-mwork*Rgas* (dVcs/Tk+dVe/Th) / (Vcs/Tk+Vk/Tk+Vr/Tr+Vh/Th+Ve/Th) A2 ; 

%working  space  pressure  time  derivative  (Pa/s) 

dmcs= (pcs*dVcs+Vcs*dpcs ) / (Rgas*Tk ) ;  %compression  space  mass  flow  (kg/s) 

dmck=-dmcs;  %compression-cooler  boundary  mass  flow  (kg/s) 

dmk=Vk*dpcs/ (Rgas*Tk )  ;  %cooler  mass  flow  (kg/s) 

dmkr=dmck-dmk;  %cooler-regen .  boundary  mass  flow  (kg/s) 

dme= (pe*dVe+Ve*dpcs ) / (Rgas*Th) ;  %expansion  space  mass  flow  (kg/s) 

dmh=Vh*dpcs/ (Rgas*Th)  ;  %heater  mass  flow  (kg/s) 

dmhe=dme;  %expansion-heater  boundary  mass  flow  (kg/s) 

dmrh=dme+dmh;  %regen . -heater  boundary  mass  flow  (kg/s) 

dmr=dmkr-dmrh;  %regen.  mass  flow  (kg/s) 

mcs=pcs*Vcs/ (Rgas*Tk ) ;  %compression  space  gas  mass  (kg) 
mk=pcs*Vk/ (Rgas*Tk ) ;  %cooler  gas  mass  (kg) 
mr=pcs*Vr/ (Rgas*Tr ) ;  %regen.  gas  mass  (kg) 
mh=pcs*Vh/ (Rgas*Th) ;  %heater  gas  mass  (kg) 
me=pcs*Ve/ (Rgas*Th) ;  %expansion  space  gas  mass  (kg) 
rowcs=mcs/Vcs ;  %compression  space  gas  density  (kg/mA3) 
rowk=mk/Vk;  %cooler  gas  density  (kg/mA3) 

rowr=mr/Vr;  %regen.  gas  density  (kg/mA3) 

rowh=mh/Vh;  %heater  gas  density  (kg/mA3) 

rowe=me/Ve;  %expansion  space  gas  density  (kg/mA3) 

uk=dmk/ (rowk*ACHx) ;  %cooler  gas  velocity  (m/s) 
ur=dmr/ (rowr*Ar ) ;  %regen.  gas  velocity  (m/s) 
uh=dmh/ (rowh*AHHx) ;  %heater  gas  velocity  (m/s) 

Rek=rowk*abs (uk) *CHxDh/muk;  %cooler  gas  Reynolds  number 
if  Rek<=1.502113290507442e+003 

frk=64/Rek;  %cooler  gas  friction  factor 

else 

frk=0 . 184/RekA ( 0 . 2 ) ;  %cooler  gas  friction  factor 

end 

Fk=-Vk* ( f rk/CHxDh+Kk/CHxLg) *rowk*uk*abs (uk) /2 ;  %cooler  gas  friction  (N) 
dpk=Fk/ACHx;  %cooler  gas  pressure  drop  (Pa) 

Dissk=abs (dpk*dmk/rowk ) ;  %cooler  gas  dissipation  (W) 

Pek=Rek*Prk;  %cooler  gas  Peclet  number 
if  Rek  <  2300 

if  Pek  <  1.5 

Nuk=4 . 187* (1-0 . 0439*Pek) ;  %cooler  gas  Nusselt  number 

else 

Nuk=3 . 6568* ( 1+1 . 22 7/ (PekA2 ) ) ;  %cooler  gas  Nusselt  number 

end 

else 

Nuk=0 . 036*RekA0 . 8* (CHxLg/CHxDh) A ( - . 055 ) *PrkA0 . 33 ;  %cooler  gas  Nusselt 
number 
end 

Rer=rowr*abs (ur ) *rDh/mur ;  %regen.  gas  Reynolds  number 
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frr=192/Rer+4 . 53*RerA (-0 . 067) ;  %regen.  gas  friction  factor 

Fr=-Vr * ( f rr /rDh+Kr /rLg) *rowr *ur *abs (ur ) /2 ;  %regen.  gas  Friction  (N) 

dpr=Fr/Ar;  %regen.  Pressure  drop  (Pa) 

Dissr=abs (dpr *dmr /rowr ) ;  %regen.  dissipation  (W) 

hAr= ( 1+1 . 16* (Rer *Prr ) A0 . 66*por A2 . 61 ) *kHer *Amr/rDh;  %regen.  Heat  transfer 
coefficient  (W /  (mA2*K)  ) 

Reh=rowh*abs (uh) *HHxDh/muh;  %heater  Reynolds  number 
if  Reh<=l ,502113290507442e+003 

frh=64/Reh;  %heater  friction  factor 

else 

frh=0 . 184/RehA ( 0 . 2 ) ;  %heater  friction  factor 

end 

Fh=-Vh* ( f rh/HHxDh+Kh/HHxLg) *rowh*uh*abs (uh) /2 ;  %heater  friction  (N) 
dph=Fh/AHHx;  %heater  pressure  drop  (Pa) 

Dissh=abs (dph*dmh/rowh) ;  %heater  dissipation  (W) 

Peh=Reh*Prh;  %heater  Peclet  number 
if  Reh  <  2300 

if  Peh  <  1.5 

Nuh=4 . 187* (1-0 . 0439*Peh) ;  %heater  Nusselt  number 

else 

Nuh=3 . 6568* ( 1+1 . 22 7/ (PehA2 ) ) ;  %heater  Nusselt  number 

end 

else 

Nuh=0 . 036*RehA0 . 8* (HHxLg/HHxDh) A (- . 055 ) *PrhA0 . 33 ;  %heater  Nusselt  number 

end 

pe=pcs+dpk+dpr+dph;  %expansion  space  pressure  (Pa) 

hAk=Nuk*kHe*Awk/CHxDh;  %heater  heat  transfer  coefficient  (W/ (mA2*K) ) 
hAh=Nuh*kHe*Awh/HHxDh;  %cooler  heat  transfer  coefficient  (W/ (mA2*K) ) 
dWcs=pcs*dVcs ;  %compression  space  PV  power  (W) 
dWe=pe*dVe;  %expansion  space  PV  power  (W) 
dW=dWcs+dWe;  %total  PV  power  (W) 

Diss=Dissh+Dissr+Dissk;  %total  dissipation  (W) 

STr= (hAr/Amr ) / (rowr*abs (ur ) *cp) ;  %regen.  Stanton  number 
NTUr=STr * (Amr /Ar ) /2 ;  %regen.  Number  of  transfer  units 
ref f =NTUr / ( 1+NTUr ) ;  %regen.  effectiveness 

dQr=Vr*cv*dpcs/Rgas+cp* (dpcs* (Vk+Vcs+Vh+Ve) +dW) /Rgas;  %regen.  Heat  transfer 
(W) 

dQk=-Vcs*dpcs ;  %cooler  heat  transfer  (W) 
if  dmk  <  0 

dQk=dQk-abs (dQr ) /dQr*abs (dQr) A (1-reff ) ;  %cooler  heat  transfer  (W) 

end 

dQh=-Ve*dpcs ;  %heater  heat  transfer  (W) 
i f  dmh  >  0 

dQh=dQh+abs (dQr ) /dQr *abs (dQr ) A ( 1-ref f ) ;  %heater  heat  transfer  (W) 


end 


